A computational framework for testing hypotheses of the minimal mechanical requirements for cell aggregation using early annual killifish embryogenesis as a model

Introduction: Deciphering the biological and physical requirements for the outset of multicellularity is limited to few experimental models. The early embryonic development of annual killifish represents an almost unique opportunity to investigate de novo cellular aggregation in a vertebrate model. As an adaptation to seasonal drought, annual killifish employs a unique developmental pattern in which embryogenesis occurs only after undifferentiated embryonic cells have completed epiboly and dispersed in low density on the egg surface. Therefore, the first stage of embryogenesis requires the congregation of embryonic cells at one pole of the egg to form a single aggregate that later gives rise to the embryo proper. This unique process presents an opportunity to dissect the self-organizing principles involved in early organization of embryonic stem cells. Indeed, the physical and biological processes required to form the aggregate of embryonic cells are currently unknown.

Methods: Here, we developed an in silico, agent-based biophysical model that allows testing how cell-specific and environmental properties could determine the aggregation dynamics of early Killifish embryogenesis. In a forward engineering approach, we then proceeded to test two hypotheses for cell aggregation (cell-autonomous and a simple taxis model) as a proof of concept of modeling feasibility. In a first approach (cell autonomous system), we considered how intrinsic biophysical properties of the cells such as motility, polarity, density, and the interplay between cell adhesion and contact inhibition of locomotion drive cell aggregation into self-organized clusters. Second, we included guidance of cell migration through a simple taxis mechanism to resemble the activity of an organizing center found in several developmental models.

Results: Our numerical simulations showed that random migration combined with low cell-cell adhesion is sufficient to maintain cells in dispersion and that aggregation can indeed arise spontaneously under a limited set of conditions, but, without environmental guidance, the dynamics and resulting structures do not recapitulate in vivo observations.

Discussion: Thus, an environmental guidance cue seems to be required for correct execution of early aggregation in early killifish development. However, the nature of this cue (e.g., chemical or mechanical) can only be determined experimentally. Our model provides a predictive tool that could be used to better characterize the process and, importantly, to design informed experimental strategies.

1 Introduction

Annual killifish have a unique early developmental pattern that differs from most teleost species. Unlike non-annual species, for whom most morphogenetic movements are concomitant, annual killifish have separated epiboly from embryo formation, resulting in an initial phase of cell dispersal that is followed by a process of cell aggregation, with embryonic cells occurring through active cell migration in the confined space between the enveloping layer (EVL) and the yolk syncytial layer (YSL) (Figure 1) (Dolfi et al., 2014; Concha and Reig, 2022). The embryonic cells remain undifferentiated, with stem-like properties during epiboly and until the end of the dispersion phase, when the cells start to aggregate at one pole of the embryo and initiate the genetic and morphogenetic processes leading to the formation of germ layers and the establishment of the embryonic axis (Wourms, 1972; Pereiro et al., 2017; Márquez et al., 2019) (Figure 1). These processes, which span between several hours and a few days depending on the specific killifish species and the environmental conditions (Dolfi et al., 2014; Márquez et al., 2019), occur in the context of a nearly spherical egg (about 1 mm diameter) and can be easily visualized in a living animal, as eggs are optically transparent and develop outside the mother (Wourms, 1972; Pereiro et al., 2017; Reig et al., 2017; Concha and Reig, 2022). The dispersion phase is characterized by a random walk of embryonic cells moving at a very low density (Márquez et al., 2019), while the cellular processes and morphogenetic mechanisms that form the aggregate are still unknown. It has been proposed that self-organizing processes may break the initial symmetry of the embryo and initiate the aggregation process (Pereiro et al., 2017; Abitua et al., 2021), since the molecular signals involved in embryo formation are apparently non-polarized during the stages prior to the aggregate formation. However, it cannot be ruled out that an organizing center, possibly located in extraembryonic structures, provides the signals that initiate the aggregate formation (Pereiro et al., 2017) as has been shown in other non-annual teleost species (Concha and Reig, 2022).

FIGURE 1

FIGURE 1. Stages of the early development of annual killifish as recorded in vivo. Cartoon representation (A) and confocal images of embryonic cell nuclei (B) stained using the FUCCI construct, illustrating the stages of the early development of killifish embryos. Following the dispersive state (left), the cells move directionally toward a pole of the embryo (center) and form an aggregate (right). The dashed red circles indicate the approximate outline of the embryo. The time between the stages varies between several hours to a few days, depending on the environmental condition and specific killifish species (Dolfi et al., 2014; Márquez et al., 2019). Scale bar, 200 μm.

In silico modeling proved to be a powerful tool to accurately capture the essential features of various biological systems, such as wound healing (Ravasio et al., 2015a), tissue expansion (Ravasio et al., 2015b), cancer invasion (Stichel et al., 2017), and embryonic development (CONTE et al., 2008; Cai et al., 2016; Pereiro et al., 2017; Stepien et al., 2019). Thus, it has been proposed that they could be used as predictive tools to design informed experimental strategies (Kabla, 2012; Phillips, 2015). Here, we used an in silico model to understand the mechanical requirements for killifish cells to 1) remain in a dispersed state and 2) aggregate at the embryo pole to initiate embryogenesis. In our model, motile cells are represented by three-dimensional self-propelled particle spheres, which is a 3D framework commonly used to model collective cell dynamics (Belmonte et al., 2008; Henkes et al., 2011; Sepúlveda et al., 2013; Méhes and Vicsek, 2014; Tarle et al., 2015). This approach has been extended to incorporate biologically relevant interactions such as cell–substrate friction, intercellular and cell–substrate adhesions (Kanchanawong et al., 2010; Bertocchi et al., 2017), and contact inhibition of locomotion (CIL) (Abercrombie and Heaysman, 1954; Roycroft and Mayor, 2016). When generalized, this model can exhibit diverse dynamic states, such as gas phases, polar liquids, and 3D aggregates, depending on the parameter explored (Bray, 1993; Mladek et al., 2006; Moreno and Likos, 2007; Redner et al., 2013). Although these states were experimentally observed at high cell densities, it is an open question as to whether such mechanisms could account for the aggregation behavior observed at low cell densities found in the early stages of killifish development. Furthermore, to date, models have not considered the specific geometry of this process, such as the cells moving in confinement and on curved surfaces with spherical topology. Typically, these studies use periodic boundary conditions on a plane, giving a toroidal topology (Bellomo et al., 2015). Our model incorporates realistic conditions in terms of cell density, geometrical and mechanical properties of the EVL and the YSL, and their effect on the dynamics of embryonic cells. Thus, the in silico investigation presented aims to provide a flexible framework to model the early teleost development, which can help predict the minimal mechanical requirements for cell aggregation under the specific conditions of annual killifish early development. As the in vivo system is poorly understood and presents various intrinsic experimental challenges (e.g., coriaceous chorion), our forward engineering approach can provide useful information, enabling an informed experimental investigation of the biological system.

2 Results

During the early stages of killifish embryo development, undifferentiated stem cells move tangentially between the inner surface of the epithelial enveloping cell layer and the yolk syncytial layer (Figure 2A). The system is modeled here in two distinct ways: a cell-autonomous system that includes mechanisms that are intrinsic to the cells and the same cells that are under the influence of guidance from the environment toward an organizing center.

FIGURE 2

FIGURE 2. Schematic representation of the model’s major components. (A) Embryonic cells move between the EVL and the YSL, which are roughly spherical. They are subject to forces due to the deformation of the adjacent cell layers. (B) Polarity of the cells is constrained to the tangential plane as they move. (C) Contact inhibition of locomotion repolarizes cells away from the average position of their neighbors, and adhesion forces attract them to their neighbors. (D) Angle with respect to the average neighbor position. The dashed vectors indicate the directions to the neighbors of the red cell, the bottom solid arrow indicates the reference direction, and ${\theta }_i^f$ is the angle to which the cells repolarize due to CIL. The direction of the adhesion force is shown for reference.

2.1 Cell-autonomous system

We present, here, a variation of the model proposed by Smeets et al. (2016) (Basan et al., 2013) for autonomous motile cells, extended to three dimensions and including the physical and geometrical constraints imposed by the EVL and the YSL, where the cells migrate tangentially to the surface of the YSL (Figure 2B).

2.1.1 Equations of motion

The cells have an intrinsic motile force, $F_m$, that drives them forward along their direction of polarity ${\widehat{\mathit{p}}}_{\mathit{i}}$ (Figure 2C). The cells are subject to viscous forces from the substrate with the coefficient ${\gamma }_s$ and from other cells $F_{ij}^{cc}$. The equation of motion is as follows:

$F_m{\widehat{\mathit{p}}}_i={\gamma }_s{\dot{\mathit{x}}}_i+{\displaystyle \sum }_j{F}_{ij}^{cc}{\widehat{\mathit{n}}}_{ij}(1)$

with the left-hand term and the motile force with direction ${\widehat{\mathit{p}}}_i$, where ${\gamma }_s$ is the substrate viscosity and $F_{ij}^{cc}$ is the force between the cells $i$ and $j$, acting in the normal direction to their surface ${\widehat{\mathit{n}}}_{\mathit{i}\mathit{j}}$ at the point of contact. The normal direction at the point of contact in the case of the two spheres is simply the direction between their two centers so that:

${\widehat{\mathit{n}}}_{\mathit{i}\mathit{j}}=\frac{\left({{\mathit{x}}_{\mathit{i}}-\mathit{x}}_j\right)}{d_{ij}},(2)$

with $d_{ij}=\left|{{\mathit{x}}_i-\mathit{x}}_j\right|$ being the distance between cells. The force $F_{ij}^{cc}$ is defined in terms of $d_{ij}$.

2.1.2 Cell–cell adhesion and repulsion

The adhesion and repulsion between cells may play a significant role in the formation of aggregates. These forces are determined by the interaction between the cell–cell adhesion energy $W_c$ and the cell–substrate adhesion energy $W_s$, such that the intercellular force is given by (Smeets et al., 2016):

$F_{ij}^{cc}=\frac{2}{R}\left[W_s-\frac{W_s+W_c}{R}\left(d_{ij}-R\right)\right],(3)$

for pairs of cells i and j in contact ($d_{ij}<2R$) with the cell radius $R$. Here, $F_{ij}^{cc}=0$ when $d_{ij}\ge 2R$, as the cells are not in contact. We model the EVL and the YSL as single large spherical cells of radius $R_E$ centered on the origin:

$F_i^{YSL}=\frac{2}{R}\left[W_s-\frac{W_s+W_c}{R}\left(\left|{\mathit{x}}_{\mathit{i}}\right|-R_E-R\right)\right],$

which is the force applied to cell i by the YSL when $\left|x_i\right|-R_E<R$ ($F_i^{YSL}=0$ otherwise), and

$F_i^{EVL}=\frac{2}{R}\left[W_s-\frac{W_s+W_c}{R}\left(R_E-\left|{\mathit{x}}_i\right|-R\right)\right],$

which is the force applied to cell i by the EVL when $R_E-\left|{\mathit{x}}_{\mathit{i}}\right|<R$ ($F_i^{EVL}=0$ otherwise). In this way, the cells experience forces that tend to maintain them on a sphere of radius $R_E$. It should be noted that unlike in the work of Smeets et al. (2016), there is no cut-off of intercellular forces when the cells are closer than R, and no cells are removed from the simulation to model multiple layers. This is not necessary since our model is fully three-dimensional, and due to the spherical geometry, cells are naturally forced outward to form multilayered aggregates, which can be observed under large adhesive energies. Since the dynamics of the system depend upon the relation between cell–cell and cell–substrate energies rather than their absolute strengths, from here on, we use $W_s=1$ and simply vary $W_c$.

2.1.3 Repolarization and rotational diffusion

In our model, the cell polarity vector ${\widehat{\mathit{p}}}_{\mathit{i}}$ is constrained to the tangential plane between the EVL and the YSL, and the motile force is as follows:

${\mathit{F}}_m=F_m{\widehat{\mathit{p}}}_{\mathit{i}},(4)$

where ${\widehat{\mathit{p}}}_{\mathit{i}}={\left(\mathit{cos}{\theta }_i,\mathit{sin}{\theta }_i\right)}^T$, and $F_m$ is the magnitude of the motile force. We used the model by Smeets et al. (2016) for the repolarization of the direction vector of the cell ${\widehat{\mathit{p}}}_{\mathit{i}}$, which has angle ${\theta }_i$ in the tangential plane (Figure 2C):

${\dot{\theta }}_i=-f_{pol}\left({\theta }_i-{\theta }_i^{*}\right)+\xi \sqrt{2D_r},(5)$

where ${\theta }_i^{*}$ is the target or desired direction of the cell, $f_{pol}$ is the rate of repolarization, $D_r$ is the rate of angular diffusion, and $\xi$ is a Gaussian noise process. We may define any formulation for the desired direction depending on which type of repolarization process we wish to consider. This equation can be normalized by giving a single dimensionless parameter, $\psi =f_{pol}/2D_r$. In the following, we fix $D_r$ at unity.

2.1.4 Contact inhibition of locomotion

Contact inhibition of locomotion (CIL), which works as a repulsion interaction causing cells to steer away from each other, is another important determinant of cell migration, which we wish to test in our model. For CIL, we model the repolarization direction, ${\theta }_i^{*}={\theta }_i^f$, as the direction pointing away from the average position of each cell’s neighbors (contacting cells—Figure 2D), so if the average position of the neighboring cells is

${\overline{\mathit{x}}}_{\mathit{i}}={\displaystyle \sum }_j{\mathit{x}}_{\mathit{j}},$

for the neighboring cells j, the desired direction vector in the tangent plane with normal ${\widehat{\mathit{n}}}_{\mathit{i}}$ is

${\mathit{v}}_i^f=\left({\overline{\mathit{x}}}_i-{\mathit{x}}_{\mathit{i}}\right)-\left[\left({\overline{\mathit{x}}}_{\mathit{i}}-{\mathit{x}}_{\mathit{i}}\right)\cdot {\widehat{\mathit{n}}}_i\right]{\widehat{\mathit{n}}}_i.$

Thus, the cells will turn away from contacting cells at a rate of $f_{pol}=f_{cil}$, so that

${\dot{\theta }}_i=-f_{cil}\times \mathit{arccos}\left(\frac{{\mathit{v}}_i^f\bullet {\widehat{\mathit{p}}}_i}{\left|{\mathit{v}}_i^f\right|}\right)+\sqrt{2D_r}\xi .(6)$

The system can again be normalized by the dimensionless parameter $\psi =f_{cil}/2D_r$ since $f_{pol}=f_{cil}$.

2.1.5 Numerical simulations

This model can be solved numerically (see Methods) for a range of parameters, providing a powerful tool to simulate developmental processes occurring in a realistic geometrical setting. Using the model described previously, we first examined the effect of the cell density on the collective behavior of cells, with other parameters for cell–cell adhesion and contact inhibition being fixed at W_c = 1 and Ψ = 0.5, respectively. Spherical cells of radius 1 were initialized at random positions on the surface of a sphere with a radius 25, representing the YSL, and their migration simulated 500 time steps. At a cell packing density ($\mathrm{\Phi }=n{R}^2/4R_E^2$ for n cells on an embryo of radius $R_E$) similar to the annual killifish embryo (see Supplementary Material) (Φ = 0.2; 500 cells; Figure 3A and Supplementary Movie S1), the cells were able to aggregate into separated clusters and, as shown in previous studies for two-dimensional systems (Basan et al., 2013; Smeets et al., 2016), form a cohesive single aggregate at a high density (Φ = 0.56; 500 cells; in a smaller sphere of radius 15; Figure 3B and Supplementary Movie S2).

FIGURE 3

FIGURE 3. Aggregation in the cell-autonomous system. (A) Cell density $\varphi =0.2$ (500 cells on a sphere of radius 25), $\psi =0.5$, and $W_c=1.$ (B) Cell density $\varphi =0.56$ (500 cells on a sphere of radius 15), $\psi =0.5$, and $W_c=1$. The images on the left show the initial conditions, and the images on the right show the simulations after 500 time steps. The YSL is drawn as a transparent gray sphere; hence, the cells on the far side appear darker. The three-dimensional perspective means that these cells also appear smaller.

To understand the effect of cell–cell adhesion and CIL on these dynamics, we scanned the parameter space, simulating the system with a range of values of W_c and Ψ for each density (Φ = 0.2 or 0.56). Figure 4 shows the final configurations of cells after 500 time steps for each parameter combination. A clear pattern emerges, where at low values of W_c, the system maintains its dispersed condition, and at high W_c, aggregation occurs (Figures 4A, B). Similar to what was shown previously, large single aggregates can be obtained at high densities and strong cell–cell adhesions, whereas at a first approximation, Ψ appears to have only a marginal effect on the qualitative appearance of the cell aggregate. To quantify these effects, we computed the maximum cluster size at the end of simulations for each of the parameter combinations (Figure 5), and we also analyzed its time evolution (Figure 6). This analysis clearly confirms that larger aggregates are formed at a higher W_c, where both cell densities are considered, and that W_c has a greater influence on the formation of aggregates as compared to Ψ. Interestingly, it also showed a sort of weak biphasic effect of Ψ on the cluster size, with larger aggregates being formed at intermediate values (between 0.75 and 1.5). However, high variations in final cluster sizes can be seen at these values, which may be due to the more pronounced stochasticity of the process.

FIGURE 4

FIGURE 4. Parameter scan of the cell-autonomous system over different cell adhesion and CIL parameters. (A) Cell density of $\varphi =$ 0.2 (500 cells on a sphere of radius 25). (B) Cell density of $\varphi =$ 0.56 (500 cells on a sphere of radius 15). The dispersed states are maintained at low cell–cell adhesion, and clustering is increased by high cell–cell adhesion. High levels of CIL (Ψ) also increase cell clustering at moderate levels of cell–cell adhesion. The YSL is drawn as a transparent gray sphere; hence, the cells on the far side appear darker. The three-dimensional perspective means that these cells also appear smaller.

FIGURE 5

FIGURE 5. Maximum cluster size as a function of CIL and cell–cell adhesion. The maximum cluster sizes for the final states (after 500 time steps, each grid point represents the average of 25 simulations) of simulations with cell densities of (A) $\varphi =$ 0.2 (500 cells on a sphere of radius 25) and (B) $\varphi =$ 0.56 (500 cells on a sphere of radius 15). Clustering is increased by high cell–cell adhesion (W_c). Both high and low levels of CIL (Ψ) decrease cell clustering. The plots are the averages of 25 independent simulations.

FIGURE 6

FIGURE 6. Time dynamics of the aggregation of the cell-autonomous system measured by the maximum cluster size. The cluster size is plotted here as a function of time for different cell densities, both conditions include 500 cells in each simulation ($\varphi =0.2$, (A–C); $\varphi =0.56$, (D–F) for different $W_c$ (0, 0.5, and 1 for the first, second, and third column, respectively) and for two different values of $\psi$ (0 for the blue circles and 2 for the red ones). The circles are an average from 50 different simulations, and the dashed lines represent the range within the standard deviation. All plots are generated using a logarithmic scale.

Numerical simulations also allow us detailed insights into the dynamics of the aggregation process, which can be quantified by the time variation in the maximum cluster size (Figure 6). With no cell–cell adhesion (W_c = 0), the system was largely static, with no increase in the cluster size at either density tested after 10 dimensionless time units. This is a condition that closely resembles the dynamics observed during the dispersion state. However, as W_c was increased, the maximum size of the clusters tended toward the power law dynamics. On the other hand, CIL appeared to have a marginal effect on the dynamics (rate) of aggregation for intermediate and high values of W_c (Figures 6B, C, E, F), whereas a marked effect of CIL can be seen when W_c is absent (Figures 6A, D). In these conditions, a high CIL (Ψ = 2) caused an initial decline in the size of the cluster and a generally low aggregation as compared to Ψ = 0. This phenomenon, which was more pronounced at high cell densities, most likely reflects the scattering effect provided by CIL.

2.2 Environmental guidance

The results presented previously show that purely cell-autonomous behaviors could explain the dispersion state by keeping low cell–cell adhesion, but they were not sufficient to lead to the formation of a single aggregate at one of the embryonic poles, as seen in annual killifish early embryogenesis. It is, therefore, possible that the information provided by environmental cues is needed in the form of an organizing center that causes cells to orient toward a specific position of the embryo, where possibly a site-specific increase in W_c for cells reaching the location could initiate the aggregation process. The cells might preferentially move toward the organizing center by a variety of mechanisms including chemotaxis, durotaxis, and haptotaxis (CARTER, 1967; Bellomo et al., 2015; Espina et al., 2022). A powerful feature of our model and its software implementation is that such external environmental factors can easily be included. As a proof of concept, we show, here, a simple model of repolarization that shows a bias for the cells to repolarize toward the organizing center. However, experiments support a variety of possible guidance mechanisms (Sarris and Sixt, 2015a). As a feasibility study, in the Supplementary Material, we further expand two more models (i.e., “adjustment of directional speed along gradient” and “slowing down at the source”) to demonstrate that our modeling framework provides a robust and flexible tool to model a variety of taxis models (see Supplementary Material, Figures S1, S2).

2.2.1 External taxis

If the organizing center is at position $x_{org}$, the vector pointing from cell i to the organizing center in the tangential plane is as follows:

${\mathit{v}}_i^{org}=\left({\mathit{x}}_{org}-{\mathit{x}}_i\right)-\left[\left({\mathit{x}}_{org}-{\mathit{x}}_i\right)\cdot {\widehat{\mathit{n}}}_i\right]{\widehat{\mathit{n}}}_i.(7)$

Then, the angular equation of motion becomes

${\dot{\theta }}_i=-f_{cil}\times \mathit{arccos}\left(\frac{{\mathit{v}}_i^f\bullet {\widehat{\mathit{p}}}_i}{\left|{\mathit{v}}_i^f\right|}\right)-f_{tax}\times \mathit{arccos}\left(\frac{{\mathit{v}}_i^{org}\bullet {\widehat{\mathit{p}}}_i}{\left|{\mathit{v}}_i^{org}\right|}\right)+\sqrt{2D_r}\xi (8)$

In this simple model, the rate of repolarization does not depend on the distance from the organizing center, nor are there any effects on the speed of cell motion. Various mechanisms and models of taxis have been proposed (Sarris and Sixt, 2015b), which while not considered here, are straightforward to implement in our modeling framework.

2.2.2 Numerical simulations

To test the effect of this simple external guidance (taxis), we then performed a series of simulations using the more realistic (i.e., closest to in vivo) conditions with cell number = 500 and sphere radius = 25, for density $\varphi =0.2$, while testing the effect of varying the external taxis repolarization rate (f_tax), cell adhesion ($W_c$), and CIL ($\psi$). As expected, at a low f_tax (0.01), cells had an overall tendency to move toward the organizing center but did not form a single aggregate within the time of our simulation as the properties of the cell-autonomous system, such as rotational diffusion and CIL, prevailed (Figures 7A, D). At this low f_tax regime, the cells mostly remained dispersed on the surface of the sphere as single cells for $W_c$ = 0 (Figure 7A) or formed small aggregates that slowly coalesced at the organizing center for $W_c$ = 1 (Figure 7D). On the other hand, intermediate and high strengths of f_tax (0.1 and 1) showed features similar to the dynamics of aggregation principally depending on W_c and CIL, while f_tax determined the speed (rate) and the degree of the aggregation process with intermediate values of f_tax (Figures 7B, E) being, at a first approximation and to different degrees, a slower and attenuated version of the dynamics seen for the highest f_tax value (Figures 7C, F, 8; Supplementary Movie S3–S6). Our simulations using f_tax = 1 showed four distinct and equally interesting phenotypes. Strong adhesion ($W_c=1$) combined with strong CIL ($\psi =2$) led the cells to the formation of small aggregates that fluctuated on the surface of the sphere and eventually coalesced into a large aggregate, only at long time scales (Figure 7F, red curve; Figure 8A; Supplementary Movie S3). This was also the case for the same conditions but using f_tax = 0.1, where the only noticeable difference was the rate at which small aggregates moved toward the organizing center (Figure 7E, red curve). In similar conditions of CIL ($\psi =2$), but in the absence of adhesion ($W_c=0$), the cells moved toward the organizing center and remained in a form of dispersed single-cell dynamics within a small area close to the organizing center. However, due to high CIL, the cells did not form a tight aggregate, and the maximum cluster size plateaued at about 100 cells (Figure 7C, red curve; Figure 8B; Supplementary Movie S4). This characteristic was also seen with f_tax = 0.1, where the cells randomly collided with the temporarily small clusters. However, due to the low f_tax, the cells were dispersed over a larger area as compared to f_tax = 1, thus diminishing the probability of collision (maximum cluster size stabilized at 3–5 cells). Finally, the aggregation into a single cluster at the organizing center was achieved when CIL was not present ($\psi =0$), irrespective of the value of W_c and f_tax (Figures 7B, C, E, F, blue curves; Figures 8C, D; Supplementary Movie S5, S6). However, W_c proved to play an important role in defining the mode of migration toward the organizing center, with cells at high W_c moving as small aggregates (Figure 7F, blue curve; Figure 8C; Supplementary Movie S5), whereas, in the absence of adhesion (W_c = 0), the cells moved independently as single cells (Figure 7C, blue curve; Figure 8D; Supplementary Movie S6). Furthermore, at intermediate f_tax and in the absence of cell–cell adhesion (Figure 7B, blue circles), the cells at the border of the aggregate lacked the necessary pulling force from the organizing center to remain cohesively adherent to each other; thus, the aggregate approached but never reached the maximum cluster size of 500 cells. These last two conditions (Figures 7B, C, blue circles; Figure 8D; Supplementary Movie S6) seem to recapitulate the best in vivo observations previously reported, for example by Dolfi et al. (2014), where the cells migrate as single cells to form a loose cluster (Figure 1). However, the precise details of the migration and aggregation mechanism in annual killifish are still largely elusive. Thus, the experimental observations [as in (Dolfi et al., 2014)] need to be corroborated, expanded, and carefully analyzed before being implemented in our model.

FIGURE 7

FIGURE 7. Time dynamics of aggregation for the environmental guidance model measured by the maximum cluster size. The cluster size is plotted here as a function of time for different $W_c$ [(0 for (A–C) and 1 for (D–F)], different f_tax (0.01, 0.1, and 1 for the first, second, and third column, respectively), and for two different values of $\psi$ (0 for the blue circles, and 2 for the red ones). The data points are the average from 50 different simulations, and the dashed lines represent the standard deviation. The conditions used accurately represent the in vivo density of the cells ($\varphi =0.2$, 500 cells on a sphere of radius 25). All plots are generated using a logarithmic scale. The simulations were stopped when they reached a single aggregate (maximum cluster size of exactly 500 cells). Hence, the one cluster condition is only reached in those simulations that are interrupted at an earlier time [i.e., $\psi$ = 0 in (C, E); $\psi$ = 2 in (F), whereas for $\psi$ = 0 in (B), the curve gets very close but never reaches the single cluster condition].

FIGURE 8

FIGURE 8. Frames of different time steps for an environmental guidance system (external taxis). For all simulations, f_tax is 1 and $\varphi =0.2$ (500 cells on a sphere of radius 25). We varied $\psi$ and $W_c$ to explore the emerging aggregation dynamics and phenotypes. (A) $\psi =2$; $W_c=1.$ (B) $\psi =2$; $W_c=0.$ (C) $\psi =0$; $W_c=1.$ (D) $\psi =0$; $W_c=0.$ The YSL is drawn as a transparent gray sphere; hence, the cells on the far side appear darker. The three-dimensional perspective means that these cells also appear smaller. The images in the left, middle, and right columns represent the times 0, 50, and 500, respectively, except for C, where the final configuration is reached at time 72.5. The arrow indicates the position of the organizing center.

3 Discussion

We have developed a computational model of cell migration and aggregation mechanisms for cells trapped between two concentric spheres, as is the case of the in vivo conditions of the annual killifish early embryo development. We demonstrate that numerical simulations of this model are an essential tool to test the role of physical mechanisms based on both the cell-autonomous and environmental guidance factors in driving the complex self-organizing processes during morphogenesis. Our computational model includes the full geometry of the embryo, including the spherical topology due to the constraint between the EVL and the YSL. This model was solved numerically to build phase space maps that reveal the effect of key mechanical parameters on the aggregation process.

For cell-autonomous mechanisms, the simulated aggregate phenotype was comparable to the experimental observation of the dispersed state, which could be achieved simply if the cells had low adhesion but failed to recapitulate the aggregation into a single cluster (Figures 3–6). At high adhesion, smaller clusters fused into a larger cluster, but the type of coalescence observed in vivo was never observed under these simple conditions. Thus, we explored the possibility that in addition to cell-autonomous factors, external cues possibly arising from a putative organizing center are also at play. The environmental taxis cues that guide cell migration could, in principle, be of many kinds, possibly chemical (i.e., diffusible or substrate-bound gradients of a chemoattractant) or physical (i.e., stiffness gradients). Using our model, we simulated a simple external cue guiding the direction of the migration of cells and showed that, under some conditions, it does indeed form single localized clusters (Figures 7, 8). In particular, our simulations determined that in order to achieve the dynamic and resulting phenotype of aggregation observed in early killifish development, an organizing center and low levels of CIL are essential. In addition, it appears that adhesive forces must remain low throughout the process and possibly increase only after the aggregate has formed.

Importantly, our strategy aims at building knowledge using a bottom-up approach by making predictions starting from the minimal condition for aggregation. Thus, at present, our model is kept to the most simplistic conditions possible in order to avoid overfitting of the system. However, more detailed models of cell motility or taxis, for example, may easily be constructed using our framework, given the relevant mechanistic details and more complete experimental evidence. Another aspect that might be considered is lateral inhibition, a mechanism that implicates the cell inhibition of adjacent cell activity. We are currently exploring all these mechanisms to simulate an aggregation process that is most like the experimental phenomena seen in annual killifish. Most importantly, despite its simplicity, our model is able to recapitulate the major features of cell aggregation in early annual killifish embryo development, a nearly unique model of cell aggregation in vertebrates.

4 Methods

4.1 Numerical simulations

The equations of motion (1) were solved using an implicit Euler method, with

${\mathit{x}}_i(t+\Delta t)-{\mathit{x}}_i(t)=\frac{\Delta t}{{\gamma }_s}(-F_m{\widehat{\mathit{p}}}_{\mathit{i}}+{\displaystyle \sum }_{\mathit{j}}{F}_{ij}^{cc}\frac{\left({\mathit{x}}_i\left(t+\Delta t\right)-{\mathit{x}}_j\left(t+\Delta t\right)\right)}{\left|{\mathit{x}}_i\left(t+\Delta t\right)-{\mathit{x}}_j\left(t+\Delta t\right)\right|})(9)$

so that we solve for $x_i\left(t+\Delta t\right)$ implicitly from Eq. 9. The rotational equations of motion given by Eqs 6, 8 were solved using the Euler–Maruyama scheme, such that

${\theta }_i\left(t+\Delta t\right)={\theta }_i\left(t\right)+\Delta t\left[f_{pol}\left({\theta }_i^{*}\left(t\right)-{\theta }_i\left(t\right)\right)\right]+\sqrt{2D_r}{\xi }_i,(10)$

where ${\xi }_i$ is the normally distributed zero mean random noise with variance $\Delta t$. In all simulations, $\Delta t$ was chosen as 0.1 for cell-autonomous and external taxis models.

For each set of parameters, 50 replicates were made for cell-autonomous and 50 for external taxis models.

4.2 Software implementation

All codes and models were implemented using the open-source CellModeller multicellular modeling framework (Rudge et al., 2012).

4.3 Embryo imaging

The in vivo experiments depicted in Figure 1 were performed under the license for animal housing, breeding, and manipulation that was issued by Umwelt-und Verbraucherschutzamt der Stadt Köln, with the authorization no. 576.1.36.6.G28/13 Be. All the fish used were raised in 35-L tanks at 24°C–26°C and belonged to the N. furzeri FUCCI transgenic strain published by Dolfi et al. (2019). They were fed two to three times a day with frozen Chironomus larvae or living nauplii of Artemia salina, depending on their size. The breeders were kept in 8-L tanks with one or more boxes (9 cm × 9 cm × 4 cm) filled with 2 cm of river sand and were left to spawn eggs for 1 h. The embryos were collected by sieving the sand with a plastic net and were then embedded in 2% low melting agar. The fluorescence images shown in Figure 1 were acquired using a Leica TCS SP5-X confocal microscope and the red emission channel with 543-nm lasers. A total of 40 to 60 images per embryo were acquired at a depth distance of 6 μm. The maximum intensity projections shown have been generated using ImageJ. The contrast in the images has been adjusted, but not altered, to optimize visualization.

Data availability statement

The raw data supporting the conclusion of this article will be made available by the authors upon request.

Ethics statement

The animal study was reviewed and approved by Umwelt-und Verbraucherschutzamt der Stadt Köln. Authorization No. 576.1.36.6.G28/13 Be.

Author contributions

TR wrote the simulation software; IM-R, GY, MA-L, OG-C, and ES performed all the numerical simulations; GY, AR, and TR designed and performed the analyses; IM-R and NOR contributed with theoretical formulations; LD and AC performed and supervised the animal experiments; MC, MLC, and CB contributed with critical theoretical and analytical knowledge; AR and TR designed, directed, and supervised the project; GY, IM-R, MLC, AR, and TR wrote the manuscript. All authors contributed to commenting and writing the manuscript.

Funding

This work was funded by ANID SCIA/ACT192015, ANID FONDECYT Regular 1211598, 1210872, 1221696, 1230919, and 1190806; ANID FONDEQUIP EMQ210101, EQM210020, and EQM130051; FONDAP 15150012; ICN09_015; NCN19_170; REDES170212, and PUENTE-2022-13 from Pontificia Universidad Católica de Chile.

Acknowledgments

AR and TR are grateful to PUC/VRI, PUC IIBM, and the graduate program at IIBM for the seed funding and support. NOR is thankful for the support received from PUC IIBM.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2023.959611/full#supplementary-material

References

Abercrombie, M., and Heaysman, J. E. M. (1954). Observations on the social behaviour of cells in tissue culture. II. Monolayering of fibroblasts. Exp. Cell Res. 6, 293–306. doi:10.1016/0014-4827(54)90176-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Abitua, P. B., Aksel, D. C., and Schier, A. F. (2021). Axis formation in annual killifish: Nodal coordinates morphogenesis in absence of Huluwa prepatterning. bioRxiv. doi:10.1101/2021.04.16.440199

CrossRef Full Text | Google Scholar

Basan, M., Elgeti, J., Hannezo, E., Rappel, W. J., and Levine, H. (2013). Alignment of cellular motility forces with tissue flow as a mechanism for efficient wound healing. Proc. Natl. Acad. Sci. U. S. A. 110, 2452–2459. doi:10.1073/pnas.1219937110

PubMed Abstract | CrossRef Full Text | Google Scholar

Bellomo, N., Bellouquid, A., Tao, Y., and Winkler, M. (2015). Toward a mathematical theory of Keller–Segel models of pattern formation in biological tissues. Math. Models Methods Appl. Sci.25, 1663-1763. doi:10.1142/S021820251550044X

CrossRef Full Text | Google Scholar

Belmonte, J. M., Thomas, G. L., Brunnet, L. G., de Almeida, R. M. C., and Chaté, H. (2008). Self-propelled particle model for cell-sorting phenomena. Phys. Rev. Lett. 100, 248702. doi:10.1103/PhysRevLett.100.248702

PubMed Abstract | CrossRef Full Text | Google Scholar

Bertocchi, C., Wang, Y., Ravasio, A., Hara, Y., Wu, Y., Sailov, T., et al. (2017). Nanoscale architecture of cadherin-based cell adhesions. Nat. Cell Biol. 19, 28–37. doi:10.1038/ncb3456

PubMed Abstract | CrossRef Full Text | Google Scholar

Bray, A. J. (1993). Theory of phase ordering kinetics. Phys. A Stat. Mech. its Appl. 194, 41–52. doi:10.1016/0378-4371(93)90338-5

CrossRef Full Text | Google Scholar

Cai, D., Dai, W., Prasad, M., Luo, J., Gov, N. S., and Montell, D. J. (2016). Modeling and analysis of collective cell migration in an in vivo three-dimensional environment. Proc. Natl. Acad. Sci. 113. doi:10.1073/pnas.1522656113

CrossRef Full Text | Google Scholar

Carter, S. B. (1967). Haptotaxis and the mechanism of cell motility. Nature 213, 256–260. doi:10.1038/213256a0

PubMed Abstract | CrossRef Full Text | Google Scholar

Concha, M. L., and Reig, G. (2022). Origin, form and function of extraembryonic structures in teleost fishes. Philos. Trans. R. Soc. Lond B Biol. Sci. 377, 20210264. doi:10.1098/rstb.2021.0264

PubMed Abstract | CrossRef Full Text | Google Scholar

Conte, V., Munoz, J., and Miodownik, M. (2008). A 3D finite element model of ventral furrow invagination in the Drosophila melanogaster embryo. J. Mech. Behav. Biomed. Mater. 1, 188–198. doi:10.1016/j.jmbbm.2007.10.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Dolfi, L., Ripa, R., Antebi, A., Valenzano, D. R., and Cellerino, A. (2019). Cell cycle dynamics during diapause entry and exit in an annual killifish revealed by FUCCI technology. EvoDevo 10, 29. doi:10.1186/s13227-019-0142-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Dolfi, L., Ripa, R., and Cellerino, A. (2014). Transition to annual life history coincides with reduction in cell cycle speed during early cleavage in three independent clades of annual killifish. EvoDevo 5, 32. doi:10.1186/2041-9139-5-32

PubMed Abstract | CrossRef Full Text | Google Scholar

Espina, J. A., Marchant, C. L., and Barriga, E. H. (2022). Durotaxis: The mechanical control of directed cell migration. FEBS J. 289, 2736–2754. doi:10.1111/febs.15862

PubMed Abstract | CrossRef Full Text | Google Scholar

Henkes, S., Fily, Y., and Marchetti, M. C. (2011). Active jamming: Self-propelled soft particles at high density. Phys. Rev. E 84, 040301. doi:10.1103/PhysRevE.84.040301

CrossRef Full Text | Google Scholar

Kabla, A. J. (2012). Collective cell migration: Leadership, invasion and segregation. J. R. Soc. Interface 9, 3268–3278. doi:10.1098/rsif.2012.0448

PubMed Abstract | CrossRef Full Text | Google Scholar

Kanchanawong, P., Shtengel, G., Pasapera, A. M., Ramko, E. B., Davidson, M. W., Hess, H. F., et al. (2010). Nanoscale architecture of integrin-based cell adhesions. Nature 468, 580–584. doi:10.1038/nature09621

PubMed Abstract | CrossRef Full Text | Google Scholar

Márquez, S., Reig, G., Concha, M., and Soto, R. (2019). Cell migration driven by substrate deformation gradients. Phys. Biol. 16, 066001. doi:10.1088/1478-3975/ab39c7

PubMed Abstract | CrossRef Full Text | Google Scholar

Méhes, E., and Vicsek, T. (2014). Collective motion of cells: From experiments to models. Integr. Biol. 6, 831–854. doi:10.1039/c4ib00115j

CrossRef Full Text | Google Scholar

Mladek, B. M., Gottwald, D., Kahl, G., Neumann, M., and Likos, C. N. (2006). formation of polymorphic cluster phases for a class of models of purely repulsive soft spheres. Phys. Rev. Lett. 96, 045701. doi:10.1103/PhysRevLett.96.045701

PubMed Abstract | CrossRef Full Text | Google Scholar

Moreno, A. J., and Likos, C. N. (2007). Diffusion and relaxation dynamics in cluster crystals. Phys. Rev. Lett. 99, 107801. doi:10.1103/PhysRevLett.99.107801

PubMed Abstract | CrossRef Full Text | Google Scholar

Pereiro, L., Loosli, F., Fernandez, J., Hartel, S., Wittbrodt, J., and Concha, M. L. (2017). Gastrulation in an annual killifish: Molecular and cellular events during germ layer formation in Austrolebias. Dev. Dyn. 246, 812–826. doi:10.1002/dvdy.24496

PubMed Abstract | CrossRef Full Text | Google Scholar

Phillips, R. (2015). Theory in biology: Figure 1 or figure 7? Trends Cell Biol. 25, 723–729. doi:10.1016/j.tcb.2015.10.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Ravasio, A., Cheddadi, I., Chen, T., Pereira, T., Ong, H. T., Bertocchi, C., et al. (2015). Gap geometry dictates epithelial closure efficiency. Nat. Commun. 6, 7683–7713. doi:10.1038/ncomms8683

PubMed Abstract | CrossRef Full Text | Google Scholar

Ravasio, A., Le, A. P., Saw, T. B., Tarle, V., Ong, H. T., Bertocchi, C., et al. (2015). Regulation of epithelial cell organization by tuning cell-substrate adhesion. Integr. Biol. (United Kingdom) 7, 1228–1241. doi:10.1039/c5ib00196j

CrossRef Full Text | Google Scholar

Redner, G. S., Baskaran, A., and Hagan, M. F. (2013). Reentrant phase behavior in active colloids with attraction. Phys. Rev. E 88, 012305. doi:10.1103/PhysRevE.88.012305

CrossRef Full Text | Google Scholar

Reig, G., Cerda, M., Sepulveda, N., Flores, D., Castaneda, V., Tada, M., et al. (2017). Extra-embryonic tissue spreading directs early embryo morphogenesis in killifish. Nat. Commun. 8, 15431. doi:10.1038/ncomms15431

PubMed Abstract | CrossRef Full Text | Google Scholar

Roycroft, A., and Mayor, R. (2016). Molecular basis of contact inhibition of locomotion. Cell. Mol. Life Sci. 73, 1119–1130. doi:10.1007/s00018-015-2090-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Rudge, T. J., Steiner, P. J., Phillips, A., and Haseloff, J. (2012). Computational modeling of synthetic microbial biofilms. ACS Synth. Biol. 1, 345–352. doi:10.1021/sb300031n

PubMed Abstract | CrossRef Full Text | Google Scholar

Sarris, M., and Sixt, M. (2015). Navigating in tissue mazes: Chemoattractant interpretation in complex environments. Curr. Opin. Cell Biol. 36, 93–102. doi:10.1016/j.ceb.2015.08.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Sarris, M., and Sixt, M. (2015). Navigating in tissue mazes: Chemoattractant interpretation in complex environments. Curr. Opin. Cell Biol. 36, 93–102. doi:10.1016/j.ceb.2015.08.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Sepúlveda, N., Petitjean, L., Cochet, O., Grasland-Mongrain, E., Silberzan, P., and Hakim, V. (2013). Collective cell motion in an epithelial sheet can Be quantitatively described by a stochastic interacting particle model. PLoS Comput. Biol. 9, e1002944. doi:10.1371/journal.pcbi.1002944

PubMed Abstract | CrossRef Full Text | Google Scholar

Smeets, B., Alert, R., Pešek, J., Pagonabarraga, I., Ramon, H., and Vincent, R. (2016). Emergent structures and dynamics of cell colonies by contact inhibition of locomotion. Proc. Natl. Acad. Sci. 113, 14621–14626. doi:10.1073/pnas.1521151113

PubMed Abstract | CrossRef Full Text | Google Scholar

Stepien, T. L., Lynch, H. E., Yancey, S. X., Dempsey, L., and Davidson, L. A. (2019). Using a continuum model to decipher the mechanics of embryonic tissue spreading from time-lapse image sequences: An approximate Bayesian computation approach. PLOS ONE 14, e0218021. doi:10.1371/journal.pone.0218021

PubMed Abstract | CrossRef Full Text | Google Scholar

Stichel, D., Middleton, A. M., Muller, B. F., Depner, S., Klingmuller, U., Breuhahn, K., et al. (2017). An individual-based model for collective cancer cell migration explains speed dynamics and phenotype variability in response to growth factors. npj Syst. Biol. Appl. 3, 5. doi:10.1038/s41540-017-0006-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Tarle, V., Ravasio, A., Hakim, V., and Gov, N. S. (2015). Modeling the finger instability in an expanding cell monolayer. Integr. Biol. (Camb) 7, 1218–1227. doi:10.1039/c5ib00092k

PubMed Abstract | CrossRef Full Text | Google Scholar

Wourms, J. P. (1972). Developmental biology of annual fishes. I. Stages in the normal development ofAustrofundulus myersi Dahl. J. Exp. Zoology 182, 143–167. doi:10.1002/jez.1401820202

CrossRef Full Text | Google Scholar