# Extracellular matrix flow guides emergent in-vitro morphogenesis

Abdel Rahman Abdel Fattah, Francesca Sgualidino, Lorenzo Mattolini, Myriam Roussigne, Suresh Poovathingal, Kristofer Davie, Pieter Baatsen, Katlijn Vints, Natalia Gounko, Giulio Superti-Furga, et al.

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Abdel Rahman Abdel Fattah, Francesca Sgualidino, Lorenzo Mattolini, Myriam Roussigne, Suresh Poovathingal, et al.. Extracellular matrix flow guides emergent in-vitro morphogenesis. 2023. hal-04780733

HAL Id: hal-04780733

https://hal.science/hal-04780733v1

Preprint submitted on 13 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Extracellular matrix flow guides emergent in-vitro morphogenesis

**Abdel Rahman Abdel Fattah**

CeMM Research Center for Molecular Medicine of the Austrian Academy of Sciences, Austria

**Francesca Sgualidino**

KU Leuven

**Suresh Poovathingal**

VIB Center for Brain & Disease Research

**Kristofer Davie**

VIB Center for Brain & Disease Research https://orcid.org/0000-0003-2182-1249

**Pieter Baatsen**

VIB/University of Leuven

**Katlijn Vints**

Electron Microscopy Platform of the VIB Bio Imaging Core at the KU Leuven, Belgium

**Natalia Gounko**

Electron Microscopy Platform of the VIB Bio Imaging Core at the KU Leuven, Belgium
https://orcid.org/0000-0003-0791-9477

**Giulio Superti-Furga**

CeMM Research Center for Molecular Medicine of the Austrian Academy of Sciences
https://orcid.org/0000-0002-0570-1768

**Lorenzo Mattolini**

KU Leuven

**Myriam Roussigne**

Université Toulouse 3 Paul Sabatier

**Bertrand Benazeraf**

Université Toulouse 3 Paul Sabatier

**Adrian Ranga** (adrian.ranga@kuleuven.be)

KU Leuven

**Biological Sciences - Article**

**Keywords:**

**Posted Date:** June 9th, 2023
DOI: https://doi.org/10.21203/rs.3.rs-2957065/v1

**License:** © Ⓟ This work is licensed under a Creative Commons Attribution 4.0 International License.

Read Full License

**Additional Declarations:** There is **NO** Competing Interest.
Extracellular matrix flow guides emergent in-vitro morphogenesis

Abdel Rahman Abdel Fattah<sup>1,2*</sup>, Francesca Sgualdino<sup>1</sup>, Lorenzo Mattolini<sup>1</sup>, Myriam Roussigne<sup>3</sup>, Suresh Poovathingal<sup>4</sup>, Kristofer Davie<sup>4</sup>, Pieter Baatsen<sup>5</sup>, Katlijn Vints<sup>5</sup>, Natalia V. Gunko<sup>5</sup>, Giulio Superti-Furga<sup>2,6</sup>, Bertrand Benazeraf<sup>3</sup>, Adrian Ranga<sup>1*</sup>

**Abstract**

Morphogenesis of vertebrate embryos requires extracellular matrix (ECM) in multiple developmental contexts<fnref n="1" />. Symmetry-breaking events in initially homogenous tissues trigger collective cellular dynamics which follow fluid-like behaviors, and the ECM, which is in constant engagement with morphing tissues, also undergoes fluid-like motion<fnref n="2" />⁻<fnref n="4" />. Cell and ECM deformations require coordination, however the interplay between them at cellular and tissue-scale remains largely unexplored. Here, we reveal a novel mechanism coupling morphogenetic events and epithelia-driven ECM flow. We show that exposure to ECM components triggers periodic morphogenesis of dome-shaped structures in human pluripotent stem cell monolayers, driven by directional ECM flow towards the sites of morphogenesis. We show that this flow is initiated by local symmetry-breaking events, is driven by microvilli and requires unperturbed flow conditions and cytoskeletal contractility. A theoretical model shows that a reaction-diffusion-like mechanism is responsible for organizing local morphogenesis into global tissue-wide events. The cell patterning landscape predicted by the model is experimentally recapitulated during mesoderm differentiation. These results uncover a new role for the ECM: in addition to being a site of cell adhesion, structural support and mechanical force generation<fnref n="5" />, as is currently established, we show that directional transport of ECM, termed ECM flow, is a major contributor in sustaining epithelial morphogenesis and differentiation.

<footnote n="1"><sup>1</sup> Laboratory of Bioengineering and Morphogenesis, Biomechanics Section, Department of Mechanical Engineering, KU Leuven, Leuven, Belgium</footnote>

<footnote n="2"><sup>2</sup> CeMM Research Center for Molecular Medicine of the Austrian Academy of Sciences, Vienna, Austria</footnote>

<footnote n="3"><sup>3</sup> Molecular, Cellular & Developmental Biology Unit, Center for Integrative Biology, University Paul Sabatier, CNRS, Toulouse, France</footnote>

<footnote n="4"><sup>4</sup> Center for Brain & Disease Research, VIB-KU Leuven, Leuven, Belgium</footnote>

<footnote n="5"><sup>5</sup> Electron Microscopy Platform of the VIB Bio Imaging Core at the KU Leuven, Leuven, Belgium</footnote>

<footnote n="6"><sup>6</sup> Center for Physiology and Pharmacology, Medical University of Vienna, Vienna, Austria</footnote>

<footnote>* Corresponding authors: adrian.ranga@kuleuven.be, aabdelfattah@cemm.at</footnote>
## Main

Epithelia achieve their form in the course of embryonic development through coordinated changes in cell identity, position and shape. The ECM, and specifically the motility of ECM components, is emerging as a critical mediator of these changes in various developmental contexts. At the onset of gastrulation, convective-like ECM currents have been observed at primitive streak formation⁶. When the primitive streak has been established, but before the formation of the primitive groove, glucosamine-based ECM components are trafficked between the hypoblast and epiblast⁷. The active transport of laminin and fibronectin has been reported to be necessary for subsequent axis elongation in the caudal presomitic mesoderm³. However, the mechanisms underlying tissue-ECM interactions are not well understood.

In vitro systems which combine ECM components and stem cell-derived multi-cellular constructs offer promising experimental models to advance our understanding of the interplay between ECM and morphogenesis. The ECM in such systems can play a dual role, both in imposing biophysical constraints on embedded tissues, thereby modulating their growth⁸,⁹, as well as in acting via biochemical signaling to regulate patterning¹⁰. Depending on its effective concentration, the ECM's mechanical properties, such as viscoelasticity, can be fluid-like, allowing it not only to deform in solid-like form but also flow akin to a liquid, and has been shown to support embryo-like in vitro morphogenesis. Recent examples of the exogenous effect of ECM include the emergence of in vivo-like neural tube and somite structures¹¹ and axis elongation in mouse¹² and human¹³ gastruloids upon exposure to Matrigel, a naturally derived viscoelastic ECM, neural tube folding upon neural induction of human pluripotent stem cells (hPSC) micropatterns¹⁴, as well as mesenchymal avian dermal morphogenesis via collagen fiber alignment¹⁴. It is nonetheless unclear how feedback between ECM and tissue morphogenesis occurs, particularly during critical symmetry-breaking events in early development, and whether these are conserved mechanisms underlying early morphogenetic events³,⁶.
Here, we elucidate the interactions between ECM and epithelia which are necessary to sustain morphogenesis using an hPSC-derived model of early development. We show that periodic morphogenetic events emerge in hPSC monolayers upon exposure to Matrigel and demonstrate that this phenomenon requires ECM flow in the direction of the symmetry-breaking event. A theoretical model based on these experimental observations identifies a reaction-diffusion-like mechanism underlying the relationship between cell polarity and ECM flow divergence. This model enables us to predict tissue geometry-specific cell polarity, which we validate *in vitro* through perturbation experiments. By imaging the apical cell surface using time-lapse and scanning electron microscopy, we identify that apical microvilli motion drives ECM flow. Finally, we show that just the ECM flow is sufficient to define changes in cellular identity and patterning. We reveal using single cell RNA sequencing that exposure to ECM flow leads to an early primitive streak-like phenotype, which is absent in ECM-free conditions. Moreover, we show that during mesoderm differentiation this phenotype remains confined to sites of morphogenesis that are geometrically defined by convergent streams of ECM flow. Our *in-vitro* model of early gastrulation therefore shows that extracellular matrix flows towards sites of primitive streak-like morphogenesis and that this flow is necessary to maintain early gastrulation and differentiation into mesoderm. This study unveils a completely new role for the ECM and a novel underlying mechanism: the ECM flows towards sites of symmetry-breaking, where it accumulates and serves to maintain morphogenesis and its motion in this context is not entrained or driven by cell migration within a tissue, as has been previously shown, but is instead actively transported by cells using microvilli.

## Symmetry-breaking events trigger directed ECM flow guided by cell polarity

To investigate the effect of ECM on epithelial morphogenesis, we first developed a minimal in-vitro model system where we supply Matrigel diluted in growth medium to confluent hPSC monolayers. A Matrigel concentration of 10% generated cellular dome structures with regular periodicity (156.2 ± 42.3 µm) across the monolayer (Fig.1a, Supplementary Fig.1). This morphogenetic response was
accompanied by a change in cell polarity, as determined by the patterned expression of tight junction protein ZO1 in the domes, and by the accumulation of laminin, a major component of Matrigel, into a toroidal shape localized close to the domes (Fig.1a). Live confocal microscopy of fluorescently tagged exogenous laminin captured the dynamics of ECM accumulation from an initially homogenous distribution to a toroidal volume (Supplementary Fig.2, Supplementary Movie 1). This accumulation of laminin was accompanied by an increase in laminin receptor 67LR expression at the domes (Supplementary Fig.3), suggesting an active cellular response to the emerging ECM structures, which we refer to as “toroids”.

To enhance the time-resolved visualization and analysis of ECM activity we tagged ECM components with fluorescent particles (FPs). Fluorescent ECM toroids began to emerge over 24 hours at dome locations corresponding to changes in ZO1 expressions (Fig.1b, Supplementary Movie 2) and continued for an additional 24 hours (Supplementary Fig.4). Analysis by particle image velocimetry (PIV) revealed that ECM-FP accumulation was the result of Matrigel flow vectors directed towards dome positions, which represent flow convergence points (Fig.1b, Supplementary Movie 3). Conversely, divergence points were found between domes. No lateral flow patterns were observed in the absence of cells, (Supplementary Fig.5, Supplementary Movie 3), confirming that ECM flow is a cell-driven phenomenon occurring at the cell-ECM interface.

We next investigated how ECM flow was initiated. A localized symmetry-breaking event, characterized by high ZO1 polarization in the monolayer (Fig.1c), was seen to precede any ECM flow initiation. While local fluctuations in flow velocities were evidenced, the establishment of time-averaged coherent flow allowed ECM-rich toroids to form (Supplementary Fig.6). The initially slow flow accelerated to velocities of 0.2-0.3 µm/min during a transition period between 6 and 12 h after Matrigel introduction (Fig.1d). During this period, continuous ECM filaments (Fig.1c) were observed to be pulled by and towards cells surrounding the domes, producing radial changes in the velocity field immediately around the domes (Fig.1d).
To verify whether observed bulk ECM motion was guided by changes in cell polarity, we introduced fluorescently labelled iron oxide nanoparticles into our culture system to visualize local interactions between the ECM and the cell surface. (Fig.1e). Upon Matrigel addition, fluorescent clusters (FCs) converged towards the domes (Fig.1e, Supplementary Fig.7a) following linear trajectories (Supplementary Fig.7b), while in the absence of Matrigel and without ECM flow, the trajectories of FCs were random (Fig.1e, Supplementary Fig.7b). We identified a binary switch response to Matrigel, with concentrations of >0.1% resulting in toroid formations at dome locations (Supplementary Fig.8a) with constant periodicity (Supplementary Fig.8b). All Matrigel concentrations above 1% generated similar dome height (Supplementary Fig.8c) and nanoparticle trajectory length (Supplementary Fig.9), indicating that ECM flow follows a thresholded response. Surprisingly, despite the lack of active nanoparticle accumulation at Matrigel concentration of 0.1%, dome-like structures with large periodicity were nevertheless observed (Supplementary Fig.8b), suggesting that this condition represents the lower limit of cellular interpretation of the ECM microenvironment. These results suggest that relatively small amounts of ECM, below the concentrations commonly used in current in-vitro model systems, are sufficient to trigger tissue-wide dome morphogenesis and ECM flow.

We then assessed whether ECM flow was the result of cell motion. Although initial (<6h) ECM flow could be partially attributed to cell motion, ECM velocities after 6h far exceeded cell velocities, which remained random (Supplementary Fig.10a, 10b). Moreover, FCs were observed to freely traverse over cells as they converged towards the morphing dome upon exposure to Matrigel, while they remained stationary in control cases (Fig.1f). FC motion followed a radial gradient of ZO1 expression, which was highest immediately around the domes (Fig.1g). Cell coherence, a measure of cellular alignment, gradually increased from ~6h to reach a maximum at 24 h (Fig.1g, Supplementary Fig.11, 11b). Since this time period corresponds with increased flow velocities, these results suggest that tissue-level organization contributes to flow enhancement. Overall, the emergence of local
symmetry breaks in collective cell polarity triggers and guides ECM flow, resulting in global
morphogenetic activities and patterned ECM accumulation.

ECM flow: A morphogenesis feedback mechanism

To determine the relation between ECM flow and changes in cell polarity, i.e. regions of morphogenetic activity, we calculated the spatial variation of the velocity field and found that high flow divergence correlated with low ECM accumulation and constant ZO1 monolayer-like polarity, while low divergence corresponded with dense ECM formations and changes in polarity (Fig.2a). Knowing that polarity guides ECM flow, we investigated whether flow was necessary for morphogenesis by building a theoretical model based on 2-D reaction-diffusion (Fig.2b). Regions of high polarity P (as conveyed by ZO1 expression) and low ECM flow divergence E resulted in patterns representing *in silico* analogues of *in vitro* domes (Fig.2b). To quantify the accuracy of this model, we introduced *in vitro* tissue geometries to our theoretical model and compared ZO1 expressions to their *in silico* P signal analogues. Strikingly, similarities were evidenced between these *in vitro* and *in silico* polarity signals, both in overall configuration as well as in spatial distribution (Fig.2c, Supplementary Movie 4) with *in silico* dome positions only ~75 µm on average away from their *in vitro* counterparts (Supplementary Fig.12a-d, Fig.2d). Furthermore, we show that this theoretical framework linking polarity and ECM flow divergence is widely applicable for a variety of tissue geometries, with our *in vitro* validation demonstrating this approach as a powerful predictive tool for morphogenesis (Fig.2e).

Overall, *in silico* model results suggest that the relationship between ECM-flow and polarity
involves feedback, which we further tested by performing *in silico* parameter perturbations targeting
polarity signals P and flow divergence E, followed by *in vitro* verification. Our model predicts that
abrogating flow divergence results in no flow and consequently no dome formation (Fig.2f). To
experimentally reduce flow divergence *in vitro*, we introduced shaking to the culture system, resulting
in quasi-parallel non-convergent flowlines (Fig.2f) and generating monolayer-like ZO1 expressions with
few and small symmetry breaks and minimal ECM accumulation (Fig.2f, white arrowheads, **Supplementary Fig.13**), confirming *in silico* predictions.

Similarly, the polarity field was perturbed by reducing its signal transmission (Fig.2g), which was predicted to lead to no ECM flow and unrestricted high polarity values across the tissue. As predicted by simulations, dome morphogenesis was highly sensitive to inhibition of contractility and modulation of cell polarity by blebbistatin (Supplementary Fig. 14), resulting in dome-like ZO1 expressions throughout the tissue with no FP accumulations (Fig.2g).

We further tested the importance of ECM viscoelastic properties on dome morphogenesis. When Matrigel was replaced by either polyethylene glycol (PEG) hydrogel, a linear elastic matrix which cannot flow, or by PEG-Matrigel mixture, dome morphologies and ECM accumulation were abrogated (Fig.2h). Lastly, by supplementing media with only purified laminin proteins, dome morphogenesis and toroid formations were generated, indicating that laminin alone can drive morphogenesis without requiring any of the other Matrigel constituents (Supplementary Fig. 15).

Overall, our results suggest that ECM alone is insufficient to promote morphogenetic events, which instead require strongly coupled feedback between ECM flow and morphogenesis. While ECM flow is triggered by symmetry-breaking events, the maintenance and maturation of such events rely on the establishment of the flow.

## Changing cell polarity drives ECM Flow

To examine the mechanism by which ECM components are actively transported, we imaged at high spatiotemporal resolution the cells transitioning from monolayer- to dome-cells which lie in the region identified to reel-in ECM material towards the domes. Longitudinal observations (30 min) revealed that the ECM velocity field is heterogenous (Fig.3a). High ECM velocity was associated with cells undergoing changes in polarity at the dome-monolayer interface without cell translation (Fig.3b solid and dashed circles) identifying the beginning of monolayer to dome transition. Bursts of ECM motion towards the dome were also associated with this transient process, suggesting that when the
polarity change subsides, ECM flow is delayed but is initiated elsewhere where another cell's polarity is beginning to change (Fig.3b, Supplementary Movie 5). ECM flow is therefore associated with cells that are reinterpreting their mechanical microenvironment and undergoing polarity change at the dome-monolayer interface (Fig.3c).

Since directed ECM flow occurs at the cell surface, we tested whether changes in cellular microstructures such as microvilli are involved in ECM transport. Microvilli can be triggered by flow-driven shear forces¹⁵, are linked to cytoskeleton remodeling¹⁶,¹⁷, exhibit active motion¹⁸,¹⁹ and are involved in apicobasal reorientation²⁰. We used scanning electron microscopy (SEM) to reveal sparse microvilli populations of arbitrary orientation in control conditions (Fig.3d). Upon exposure to Matrigel, dome cells were devoid of microvilli on their cell culture medium-facing surface, in stark contrast to the dense populations observed on surrounding monolayer cells (Fig.3d). At the dome-monolayer interface, microvilli were oriented towards the dome and interacted with Matrigel filaments (Fig.3d red arrowhead) through bending (Fig.3d blue arrowhead).

We further investigated the surface dynamics of cells at the dome-monolayer interface through longitudinal observations of F-actin expression, which is known to anchor to microvilli¹⁸,¹⁹ (Fig.3e, Supplementary Movie 6). Exposure to Matrigel resulted in similar microvilli velocities as in control conditions (Fig.3f), with velocity magnitudes matching those previously reported¹⁹, but had clear linear and non-stochastic trajectories (Fig.3e,g) with less directional variance (Fig.3h). Inhibiting EZRIN, a linker protein between F-actin and the microvillus structure, abrogated dome morphogenesis (Fig. 3g) and prevented FC accumulation (Fig.3k). These results suggest that stochastic cell surface processes such as microvilli motion gain coherent activity when exposed to ECM components, which serves to drag ECM to sites of morphogenesis (Fig.3i).

ECM flow-morphogenesis interplay defines the patterning landscape

To investigate whether ECM flow could mediate the specification of cell fates, we first
performed scRNAseq on hPSCs in control and 10% Matrigel conditions after 24 h of Matrigel exposure
(Fig.4a). Based on specific marker genes, we identified proliferative and non-proliferative cell populations with low and high glycolytic activity as well as a transition cluster (Fig.4b). ECM exposure resulted in highly glycolytic cells, in comparison to control condition (Fig.4c). In line with our *in vitro* observations, gene ontology (GO) enrichment analysis indicated that the addition Matrigel enhanced proliferation and cytoskeletal remodeling pathways (Fig.4d).

To gain insight into the cell fates generated from hPSCs upon ECM exposure, we compared our scRNAseq dataset to that of a human gastrula²¹ (Fig.4e). Cells in both ECM and non-ECM treated conditions most closely resembled the epiblast (EP), however, more primitive streak (PS)-like and nascent mesoderm (NM)-like cells were found in the ECM-treated conditions (Fig.4e). Notably, cells in control conditions express *BMP4*, a morphogen known to inhibit PS formation²²⁻²⁴ (Fig.4e). In contrast, *NODAL*, a hallmark of early PS²⁵⁻²⁷ was expressed in ECM-treated cells (Fig.4e). In line with transcriptomic data, nuclear pSMAD1/5/9 and cytoplasmic SMAD2/3 were found in control conditions, indicating BMP but not NODAL activity (Fig.4e). Conversely, in the ECM-treated conditions, the absence of pSMAD1/5/9 but the presence of nuclear and cytoplasmic SMAD2/3 suggests a loss of BMP and gain in NODAL activity. Additionally, NODAL activity was highest in domes, and lacking in surrounding monolayers, indicating that fate changes were specifically triggered in ECM-rich regions. NANOG, a pluripotency marker but also an early PS marker²⁸ that is upregulated in the presumptive and early primitive streak²⁹,³⁰, was more abundant in ECM-treated samples (Fig.4f). NANOG expression was highest in the domes (Fig.4f), suggesting domes acquired a PS-like phenotype. Based on these observations, we identified one dome and two monolayer cell subpopulations in ECM-treated samples (Fig.4g). Pseudotime trajectory inference analysis identified an intermediate monolayer population which gave rise to both dome-like cells and monolayer cells with neuroectoderm-like (SOX2) and anterior (OTX2)³¹ identities (Fig.4g). This was underscored by an analysis showing that domes gain NANOG identity at the expense of SOX2 expression, which is retained in the monolayer (Fig.4h). Furthermore, PS markers including those related to EMT and proliferation, as well as early mesoderm marker, *MIXL1*²⁸,³², showed clear signs of localization in dome-related regions on UMAPs (Fig.4g).
These results suggests that domes, in addition to being in a PS-like state, have posterior-like as well as
mesoderm-like characteristics in comparison to the more anterior identity of the monolayer.
Therefore, despite the absence of the primitive streak and mesoderm marker T³³⁻³⁵, our data suggests
that ECM flow may prime domes for the mesoderm fate.

To assess whether mesodermal differentiation occurs preferentially at dome positions, we
exposed ECM-treated cultures to mesodermal differentiation signaling cues (Fig.4i). ECM exposure led
to increased T expression (Fig.4j). As predicted by the transcriptomic data, most T-expressing cells
were found in domes (Fig.4k). These domes, despite being less pronounced from those formed under
pluripotent conditions that lacked T expression (Supplementary Fig.16), remained sites of polarity
change and of ECM accumulation (Supplementary Fig.17). Gains in T expression occurred at the
expense of NANOG, which decreased in the domes compared to the surrounding monolayer (Fig.4l).
This follows in vivo development where high NANOG expression is eventually excluded from the PS
prior to mesoderm differentiation²⁹. Taken together, these results suggest that feedback mechanisms
between ECM flow and morphogenesis help define the early patterning landscape during human
gastrulation, likely by tuning the susceptibility of specific regions to biochemical signaling inducing cell
fate decisions.

Overall, we show that morphogenetic symmetry-breaking events encode the conversion of
polarity signals to mechanical ECM flows (Fig.4m), and that this process is mediated by local polarity
changes and coordinated microvilli motion. In doing so, a strong feedback between morphogenesis
and ECM flow is generated, which if interrupted, leads to failure of morphogenesis. Using a theoretical
model, we show this ECM-tissue interplay follows reaction-diffusion-like principles where cell polarity
anticorrelates with flow divergence in a periodic and tissue-wide emergent phenomenon. In this way,
we recapitulate and predict observed flow and morphogenesis patterns, an approach which may be
extendable to previously reported ECM-driven and micropattern-mediated morphogenetic events,
such as neural tube morphogenesis¹⁴. Indeed, such cases of organoid morphogenesis may be specific,
geometrically constrained solutions to the generalized, geometrically unconstrained ECM flow-driven
phenomena we describe here. Our observations that ECM triggers coherent microvilli motion, and that
this motion is in turn required for directional ECM flow, suggests that actin treadmilling may be
involved in linking ECM flow and changes in cell polarity. The observed expressions of NANOG and T
during dome formation models early development processes at the onset of, and during, primitive
streak formation. An increase in glycolysis during this ECM-enhanced primitive streak-like phenotype
suggests that ECM flow may also be more broadly involved in regulating metabolic changes. While bulk
ECM motion<sup>6</sup> and ECM-tissue coupling<sup>36</sup> have been reported *in-vivo*, our *in-vitro* and *in-silico* results
reveal a previously unrecognized role for ECM transport and flow feedback to support morphogenesis
and cell differentiation and patterning in epithelia. These results establish ECM flow as a novel guiding
mechanism of epithelial morphogenesis. With the increasing use of ECMs in tissue and organoid
culture<sup>11-14</sup> this approach can help unravel the role of ECM flow in guiding morphogenesis in other
model systems.
## Acknowledgements

This work was supported by FWO postdoctoral fellowship 1217220 N (A.A.), FWO grant G087018N (A.R.) KU Leuven grants C32/17/027 (A.R.) and IDN/20/007 (A.R.), and by an Allen Distinguished Investigator Award (A.R.), a Paul G. Allen Frontiers Group advised grant of the Paul G. Allen Family Foundation.

## Author contributions

A.A. designed, conducted and analyzed all experiments. F.S. conducted experiments and assisted with dome morphology measurements. L.M. M.R. B.B. and A.R. provided input on experimental design and analysis, A.A., S.P., K.D. generated the scRNAseq data and A.A. analyzed the data. P.B., A.A., K.V. and N.V.G. prepared samples for SEM, P.B. acquired SEM images and A.A. interpreted the data. A.A., G.S.F., B.B. and A.R. edited the manuscript. A.A. and A.R. interpreted the data and wrote the manuscript.
Figures
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Concomitant morphogenesis and ECM flow dynamics. **a** Schematic representation of hPSC monolayers exposed to Matrigel and Matrigel-free pluripotent media (side view). Representative images of control and Matrigel (top and side views) conditions at 24 h showing nuclei and ZO1 expression (n = 5). Scalebar 100 µm. **b** Representative 3-D reconstructed temporal dynamics of ZO1 expression and Matrigel-tagged FPs over 24 h. FPs channel color coded for height. Black arrowheads point to representative dome formation locations at 24 h in 3-D ZO1 and FPs channels. Scalebar 100 µm. Representative 3-D reconstructed image of Matrigel toroid and y-z section of FPs between two adjacent Matrigel toroids color coded for height and accompanied 24 h average velocity map color coded for velocity magnitude where white lines and arrowhead denote flow streamlines and direction (n = 4). Scalebar 50 µm. **c** Representative images showing the temporal dynamics of a representative symmetry-breaking event (ZO1 expression, purple arrowhead) and accompanied Matrigel-tagged FPs motion over 24 h at the Matrigel-monolayer interface (n = 4). Matrigel filaments indicated by white arrowheads. Scalebars 100 µm. **d** Schematic representation of dome, surrounding monolayer and ECM flow velocity</figcaption></figure>
investigation zone (side view). Instantaneous velocity profiles for each zone over 24 h color coded for velocity zones. Curve bounds represent SEM (n = 3 for a total of 40 domes per velocity zone). Schematic representation showing the temporal dynamics of a symmetry-breaking event (Z01 expression) and associated Matrigel motion. e Schematic presentation of control and 10% Matrigel conditions with fluorescent clusters (FCs) (side view). Representative images of control and 10% Matrigel conditions (n = 3). Scalebar 100 µm. f Representative images of ZO1 expression and single FC (white arrowhead) dynamics in control and 10% Matrigel conditions over 135 min. Purple arrowheads indicate symmetry-breaking event in Matrigel condition. Colored cells highlight those traversed by the FC and are color coded by order of crossing (n = 3). Scalebar 25 µm. g Schematic representations showing ZO1 expression, FC location and ZO1 expression directional coherence (side view). Kymographs describing ZO1 expression, FC location and ZO1 average expression directional coherence color coded for magnitude (n = 3 for a total of 30 domes).
<figure><img src="image_4.png" /><figcaption>**Figure 2. Dome morphogenesis and ECM flow anticorrelation.** a Representative images of ZO1 expression and FP distribution at 24 h in 10% Matrigel conditions (n = 4) and corresponding 24 h average velocity map color coded for velocity magnitude where black lines and arrowhead denote flow streamlines and direction, derived divergence field color coded for divergence magnitude with superimposed streamlines, and divergence field with superimposed ZO1 expression (white) and FPs distribution (n = 4). Scalebars 100 µm. Normalized radial profiles over half the distance between two adjacent domes of ZO1 expression, FP distribution and flow divergence magnitudes (n = 3, for a total of 50 domes). Error bands are SEM. b Schematic representation of dome morphogenesis and ECM flow relation. Schematic profile representation of ECM flow divergence (E), tissue polarity (P) and respective mathematical expressions. Representative 2-D simulation results of scaled E and P expressions color coded for magnitude. c Process of converting *in vitro* tissue layout to a simulation template. Representative comparison between simulation and *in vitro* results depicting scaled P result (n = 4 for 143 *in vitro* and 158 *in silico* domes). Scalebars 500 µm. d Average separation length between nearest *in silico* and *in vitro* domes (from c, n = 4 with 143 comparisons). Error bars are SEM. e Scaled expression of T for various domain spaces color coded for T signal concentration with representative images of *in vitro* similar results with comparable dimensions (n = 5). Scalebars 100 µm. f Schematic representation, *in silico* analysis, and *in vitro* results of failed dome morphogenesis under shaking conditions (ECM flow perturbation) (n = 3). White arrowheads indicate small symmetry-breaking events and accompanied weak ECM accumulation. Scalebars 100</figcaption></figure>
µm. g Schematic representation, *in silico* analysis, and *in vitro* results of continuous dome morphogenesis under impaired polarity conditions resulting in ECM flow abrogation (n = 3). Scalebars 100 µm. h Schematic representation of PEG (blue lines) and PEG-Matrigel conditions. Representative image of ZO1 expression and FP distribution at 24 h after matrix exposure (n = 3). Scalebars 100 µm.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Microvilli dynamics drive ECM flow. **a** Representative images of ZO1 expression and FP distribution in 10% Matrigel conditions at 24 h (n = 3). Average velocity map color coded for velocity magnitude where black lines and arrowhead denote flow streamlines and direction overlayed with ZO1 expression (n = 3). Scalebars 25 µm. **b** Time lapse of instantaneous velocity magnitude where black arrows denote local velocity vectors (n = 3). Black circles denote ZO1 expression (polarity) regions of interests. The solid circle demonstrates a decreasing trend associated with a decreasing ECM velocity magnitude in its vicinity, and the dotted circle shows an increasing trend associated with an increase in ECM velocity magnitude in its vicinity. **c** Schematic representation of cells undergoing polarity change and driving ECM flow. **d** Scanning electron microscopy images at 24 h of cells under control and 10% Matrigel images (n = 1 with technical triplicates). Scalebars 5 µm (inset scalebars 2 µm). Microvilli orientation angle color coded for angle (-90–90°) with associated representative microvilli image superimposed with orientation angle color (n = 1 with technical triplicates). Scalebars 2 µm. **e** F-actin and</figcaption></figure>
ZO1 dynamics in control and 10% Matrigel conditions at 24 h. Surface F-actin dynamics color coded for time with white dotted lines denoting individual trajectories (n = 3). Scalebars 2 µm. f Measured microvilli (F-actin) velocities over 5 min in control and 10% Matrigel conditions (n = 3 for a total of 100 F-actin foci/microvilli). Error bars are SEM. g Measured microvilli (F-actin) trajectory start-to-finish distances over 5 min in control and 10% Matrigel conditions (from f). Statistical analysis was determined by unpaired two-sided t-test p < 0.0001. Error bars are SEM. h Measured microvilli (F-actin) average trajectory angle change over 5 min in control and 10% Matrigel conditions (from f). i Top and side view representative images of ZO1 and F-actin expressions in control and 10% Matrigel conditions with FP distribution at 24 h (n = 3). Schematic representation of microvilli dynamics in control and Matrigel conditions. j Schematic representation of microvilli dynamics in 10% Matrigel with and without EZRIN inhibition (5 µM EZRINi - NSC668394) and corresponding representative images of ZO1 and EZRIN expressions (n = 3). Scalebars 100 µm. k Corresponding FCs localization in 10% Matrigel with and without EZRINi (n = 3). Scalebars 100 µm.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Single cell RNA-seq reveals the transcriptomic landscape of hPSCs in control and 10% Matrigel conditions. **a** UMAP of the combined dataset (control and 10% Matrigel). UMAP of the combined dataset with identified clusters. **b** Expression of selected marker genes for proliferation, cytoskeletal organization, extracellular matrix remodeling, cell cycling and senescence, and neuroectoderm identity in identified clusters. Expression values are normalized and centered. **c** Dataset sample proportion and sample cluster proportions. **d** GO enrichment analysis of upregulated key processes in 10% Matrigel conditions compared to control conditions include proliferation, metabolic activity, and cytoskeletal and microtubule organization. **e** Mapping of dataset to the human gastrula dataset (E-MTAB-9388). UMAPs indicating cells identified as epiblast (E)-like, primitive streak (PS)-like, and nascent mesoderm (NM)-like. Associated sample proportion for each category. UMAP denoting</figcaption></figure>
BMP4 (PS inhibitor) and NODAL (PS inducer) expressing cells. Representative images showing SMAD2/3 and pSMAD1/5/9 immunohistochemistry in control and 10% Matrigel conditions at 24 h after Matrigel exposure (n = 3). Scalebar 100 µm, inset scalebar 25 µm.

**f** Representative images showing ZO1 and pluripotency and early PS marker, NANOG, expressions in control and 10% Matrigel conditions at 24 h (n = 3) with corresponding normalized NANOG radial profiles (n = 3 for a total of 30 domes).

**g** UMAP showing identified clusters in 10% Matrigel condition. UMAP showing two distinct inferred pseudotime trajectories color coded for pseudotime value. Expression of selected marker genes associated with pluripotency, PS, ECM interaction, neuroectoderm and mesoderm differentiation accompanied by a selection of expression UMAPs, in specific clusters.

**h** Representative image showing immunohistochemistry for NANOG (primarily in the dome) and SOX2 (primarily in the surrounding monolayers) (n = 3). Accompanied scatter plot of normalized NANOG and SOX2 expressions (n = 3 for a total of 9,744 cells).

**i** Mesoderm differentiation protocol.

**j** Integrated intensity of T expression in individual cells using consecutive confocal images (z-stack) in control and 10% Matrigel conditions (n = 3 for a total of 450 cells per condition).

**k** Representative images of mesoderm differentiation endpoint day 4 expressions of ZO1 and T (n = 3) in control and 10% Matrigel conditions with corresponding normalized T radial profiles (n = 3 for a total of 30 profiles per condition).

**l** Representative images of mesoderm differentiation end point day 4 NANOG and T in 10% Matrigel condition (n = 3) with corresponding normalized NANOG and T radial profiles (n = 3 for a total of 30 domes).

**m** Proposed model of feedback interaction between ECM flow and morphogenesis.
Supplementary figures
<figure><img src="image_7.png" /><figcaption>**Supplementary Figure 1. Periodic dome formation under Matrigel conditions.** Side and top view representative images of control and Matrigel conditions at 24 h showing nuclei and ZO1 expression (n = 5). White arrowheads denote dome locations. Scalebar 100 µm.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Supplementary Figure 2. Laminin dynamics over 24 h.** Representative 3-D reconstructed images of ZO1 and laminin temporal dynamics over 24 h (n = 3). Scalebars 100 µm.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Supplementary Figure 3. Matrigel localization affects laminin receptor expression.** Representative images of ZO1, 67LR and laminin expressions in control and 10% Matrigel conditions at 24 h (n = 3). Scalebars 100 µm.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Supplementary Figure 4.** Matrigel toroid formation at 48 h. Representative 3-D reconstructed image of ZO1 expression and Matrigel-tagged FPs at 48 h (n = 2). Top view of 3-D reconstructed image of Matrigel-tagged FPs. FPs channel color coded for height. Black arrowheads point to representative dome formation event. Scalebar 100 µm.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Supplementary Figure 5. ECM dynamics in the absence of cells.** 24 h average velocity map color coded for velocity magnitude where white lines and arrowhead denote flow streamlines and direction. The recorded velocities in control conditions can be attributed to the downward ECM flow direction due to settling during Matrigel polymerization (n = 4).</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Supplementary Figure 6. FP distribution and corresponding ECM flow field.** FP distribution over 24 h with corresponding instantaneous velocity map color coded for velocity magnitude where black arrows depict velocity vector (n = 4).</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Supplementary Figure 7. Fluorescent cluster (FC) dynamics in control and Matrigel conditions.**
Representative low magnification images of control and 10% Matrigel conditions showing FCs distribution (n = 3). Scalebar 500 µm. Representative images of ZO1 expression and FC temporal dynamics with FC trajectories (white segments) after 24 h in culture. Scalebar 100 µm.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Supplementary Figure 8.** FC localization across different Matrigel conditions and effect of Matrigel concentration on periodicity and dome height. **a** Representative low magnification images FCs distribution across various concentrations (n = 3). Scalebar 500 µm. **b** Representative images of various Matrigel conditions at 24 h showing nuclei and ZO1 expression (n = 5), example dome locations are denoted by white arrowheads. Quantification of inter-dome distance at various Matrigel conditions (n = 5 for >200 domes per conditions except for 0.1% Matrigel (31 domes), statistical analysis was determined by one-way ANOVA for multiple comparisons where the 0.1% Matrigel case had *p*-value < 0.0001 compare to all other conditions. **c** Schematic representation of dome morphogenesis and height. Quantification of dome height for various Matrigel conditions (n = 3 for a total of 20 domes per condition). Statistical analysis was determined by one-way ANOVA for multiple comparisons 0.1-2% *p* = 0.0003, 0.1-4% *p* < 0.0001, 0.1-6% *p* = 0.0003, 0.1-8% *p* = 0.0005, 0.1-10% *p* < 0.0001, 1-2% *p* = 0.0029, 1-4% *p* < 0.0001, 1-6% *p* = 0.0026, 1-8% *p* = 0.0042, 1-10% *p* < 0.0001, 2-10% *p* = 0.0128, 6-10% *p* = 0.0140, 8-10% *p* = 0.0089. Representative 3-D reconstructed images showing nuclei position color coded for height (n = 3).</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Supplementary Figure 9. Fluorescent cluster (FC) start-to-finish trajectory distance.** Quantification of average trajectory start to finish distance at various Matrigel concentrations and control conditions (n = 3 for >70 trajectories per condition). Statistical analysis was determined by one-way ANOVA for multiple comparisons control-1 to 10% *p* < 0.0001, 0.1-1% *p* = 0.0063, 0.1-2 to 10% *p* < 0.0001, 1-2 to 10% *p* < 0.0001, 2-10% *p* < 0.0001. Error bars are SEM.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Supplementary Figure 10. Velocity field comparison between cell and ECM flows.** a Average flow field velocity over time for cell and ECM flows (n = 4, for multiple view fields per replicate). b Representative images of ZO1 (cells) and ECM (Fluorescent particles – FPs) and their corresponding flow fields and heatmap denoting average velocity magnitudes in 10% Matrigel conditions (n = 4).</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Supplementary Figure 11.** Cell orientation coherence dynamics in 10% Matrigel condition. **a** Representative images of ZO1 expression, ZO1 expression orientation color coded by angle and ZO1 expression directional coherence color coded by coherence magnitude (n = 3). Scalebar 50 µm. **b** Representative images showing ZO1 expression, FC localization and coherence score dynamics over 24 h (n = 3). Scalebar 100 µm.</figcaption></figure>
<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /><figcaption>**Supplementary Figure 12.** Correlation between *in silico* and *in vitro* dome positions. **a** Pooled x and y coordinates of *in vitro* and *in silico* domes (n = 4 for 143 *in vitro* and 158 *in silico* domes). **b** Heatmap showing ranked *in vitro* and *in silico* domes and the distance between them in a shared space color coded for the separation length (from **a**). **c** Correlation between ranked *in vitro* domes and respective nearest *in silico* domes, *p*-value denotes Pearson correlation analysis *p* < 0.0001 (from **b**). **d** Difference in x and y coordinates between each *in vitro* dome and the respective nearest *in silico* one (from **b**, n = 4 for 143 *in vitro* domes).</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>**Supplementary Figure 13. 2% Matrigel conditions (without shaking).** Schematic representation of dome morphogenesis and ECM flow under stable conditions, with accompanied representative streamlines and related scaled E values (n = 3). Scaled P values with corresponding ZO1 and FP distribution at 24 h in 2% Matrigel conditions.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>**Supplementary Figure 14.** *In silico* and *in vitro* contractility perturbation comparison. Representative images of scaled P results under different values of D<sub>p</sub>. Representative images of ZO1 expression results under various blebbistatin concentrations (n = 3). Scalebar 100 µm.</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>**Supplementary Figure 15. Laminin alone is sufficient to trigger dome morphogenesis and ECM toroid generation.** Representative images showing F-actin and laminin expressions in 10% laminin conditions at 24 h after laminin exposure (n = 3). Scalebar 100 µm.</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>Supplementary Figure 16. Expression of T is absent under pluripotent medium conditions.</figcaption></figure>

Representative images showing ZO1 and T expressions in control and 10% Matrigel conditions at 48 h after Matrigel exposure (day 4) (n = 3). Scalebar 100 µm.
<figure><img src="image_26.png" /><figcaption>**Supplementary Figure 17.** ECM accumulation under 10% Matrigel and mesoderm differentiation conditions. Representative images of ZO1, 67LR and laminin expressions in 10% Matrigel conditions at 48 h after mesoderm differentiation (day 4). Scalebar 100 µm.</figcaption></figure>
## Materials and Methods

**Human PSC lines.** The human PSC lines used in this study are 1) NCRM-1 (RRID:CVCL_1E71) hPSC line from NIH Center for Regenerative Medicine (CRM), Bethesda, USA. 2) ZO1 hPSC line (Mono-allelic mEGFP-Tagged TJP1 WTC, Coriell institute for Medical Research). 3) Human iPSCs (Sigma-Aldrich IPSC0028)

**Human PSC culture.** Human PSC cultures were maintained in Matrigel-coated 6 well-plates until 60-70% confluence. To pass cells, a 4 m Dispase (Sigma-Aldrich) treatment was applied at 37°C followed by PBS washes. Colonies were gently scraped prior to gentle agitation to break down the colonies. The cell concentration was diluted 1:6 and cultured in 2 mL of Essential 8 (E8) - Flex Medium Kit (ThermoFisher Scientific) supplemented with 1% Penicillin Streptomycin (GIBCO) and Y-27632 ROCK inhibitor (ROCKi) (Hellobio) at 10 µM for 24 h. The medium was replaced with 4 mL of fresh ROCKi-free E8-Flex medium for an additional 48 h when cells reached 60-70% confluence.

**Matrigel experiments.** Human hPSCs, at a confluence of 60-70%, were washed with PBS three times, followed by the application of 1 mL of TrypLE Express (GIBCO) at 37°C for 3 min for dissociation. Once dissociated into single cells, 9 mL of DMEM/F12 medium supplemented with 20% FBS (GIBCO) was added for neutralization and TrypLE Express washing. The cells were then centrifuged at 300 RCF for 3 min. After discarding the supernatant, a second wash was applied by adding 10 mL of DMEM/F12 medium containing 20% FBS followed by centrifugation at 300 RCF for 3 min. The supernatant was discarded and the cells were resuspended in 1 mL E8-Flex medium supplemented with 10 µM of ROCKi. A cell count was performed and a volume was added containing the appropriate amount of cells to cell-free medium (E8-Flex - 10 µM ROCKi) to obtain a cell density of 150,000 cells/mL. This amounts to 30,000 cell for 200 µL of media added to a 96 well-plate. The cells were left to incubate for 24 h. Next, the medium was changed fully with fresh ROCKi-free E8-Flex medium and the cell were left to incubate
for an additional 24 h. At 48 h, with confluent wells, the cells are ready for Matrigel introduction.
Matrigel was prepared and thawed on ice. Depending on the desired Matrigel concentration, pre-
cooled (on ice) ROCKi-free E8-Flex medium was supplemented with 10%, 8%, 6%, 4%, 2%, 1% or 0.1%
Matrigel. After properly mixing the Matrigel, 200 µL of the medium was added immediately to the cells
and was kept for 24 h or 48 h. In all experiments Matrigel was introduced to the cells after 48 h of
seeding. Control conditions also receive pre-cooled ROCKi-free E8-Flex medium to minimize
confounding factors.

Differentiation experiments using Matrigel utilized mesoderm differentiation medium. The mesoderm differentiation medium was comprised of DMEM/F12 (GIBCO), 1% N2 (GIBCO), 2% B-27 (GIBCO), 1 mM sodium pyruvate MEM (GIBCO), 1 mM glutamax (GIBCO), 1 mM non-essential amino acids (GIBCO) and 2% Penicillin Streptomycin (GIBCO), 10 µM of CHIR99021 (Tocris) and 500 nM of LDN193189 (Stemgent), and 20 ng/mL of fibroblast growth factor 2 (FGF2) (Peprotech). At the Matrigel addition timepoint (48 h), the mesoderm medium was supplemented with Matrigel and then added to cells for an additional 48 h. The cells were then fixed, stained with mesoderm markers, imaged and analyzed. Control conditions did not receive Matrigel. All mesoderm experiments were paralleled with control and Matrigel conditions but using pluripotent medium, which were fixed and analyzed in the same way as the mesoderm differentiation case.

For inhibition experiments, 10 µM blebbistatin (Sigma-Aldrich), or 5 µM of EZRIN inhibitor
NSC668394 (Sigma-Aldrich) was supplemented to E8-Flex medium and introduced to cells after a full
medium change at day 1 for a 24 h pretreatment prior to Matrigel addition. After 24 h, similar inhibition
media were created, supplemented with 10% Matrigel and introduced to the cells for an additional 24
h. Control conditions did not receive Matrigel.

**Laminin labeling and live tracking.** To fluorescently tag laminin present in Matrigel, while on ice, a 10% Matrigel E8-Flex medium was created. Anti-Laminin antibody (Thermofisher PA1-16730, PA1-16730) was added to a final concentrated of 1:500 and incubated on ice for 5 min. Next donkey-anti mouse
Alexa Fluor 647 (Invitrogen) was added to the medium to a final concentration of 1:1000 and incubated on ice for 5 min. 200 µL of the medium was then introduced to the cells cultured in a well of a 96 well-plate and incubated for 5 min at 37°C and 5% CO₂. Confocal microscopy immediately followed to track laminin dynamics and ZO1 expression for 24 h.

**Fluorescent particles Matrigel tagging and live tracking.** While on ice, a 10% Matrigel E8-Flex medium was created. The medium was supplemented with fluorescent particles (FluoSpheres Carboxylate-Modified Microspheres, 200 nm, Invitrogen) to a final concentration of 1:500 v/v. The medium was then introduced to the cells and incubated for 5 min at 37°C and 5% CO₂. Confocal microscopy immediately followed to track bead displacements and ZO1 expression for 24 h or 48 h.

**Shaking experiment.** 150,000 cell were seeded into a plate of a 12 well-plate following standard Matrigel experiment protocol for media changes. At day 2, a 2% Matrigel E8-Flex medium tagged with fluorescent particles was introduced to the cells, the plate was immediately placed on an orbital shaker at 40 RPM and incubated for 5 min at 37°C and 5% CO₂ for 24 h. Control Matrigel conditions were not placed on the shaker. Matrigel free conditions were present in both conditions.

**Polyethylene glycol hydrogel Experiments.** Nondegradable polyethylene glycol (PEG) hydrogel was prepared with a 0.7 kPa stiffness as previously described⁸ and supplemented with 1:500 Fluorescent particles. Here, the cell volume fraction was replaced by milliQ water. 50 µL of the unpolymerized hydrogel was introduced to the cells, instead of Matrigel, and left for 20 m at room temperature to complete the polymerization process. After which 200 µL of media was added. To supplement the hydrogel with Matrigel, the same hydrogel premixture was used and Matrigel was added to a final concentration of 10% and replaced the same volume of MilliQ water. The combined PEG-Matrigel mixture was added to the cells and left for 20 m at room temperature to complete the polymerization process. Confocal microscopy immediately followed to track the bead displacement and ZO1
expression.

**Laminin treatment experiment.** Cells were culture as in the Matrigel experiments in 96 well-plates. At 48 h, laminin (CORNING) was thawed on ice and added to ROCKi-free E8-Flex medium to a final concentration of 0.14 mg/mL (10% v/v). The medium was immediately added to the cells which were incubated for an additional 24 h. Laminin-free and Matrigel controls were also obtained for comparison.

**Fluorescently labeled Iron oxide clusters (FCs).** 50 mg of Fe₃O₄ nanoparticles (Sigma-Aldrich, 20-40 nm) were first weighed under sterile conditions and added to 1 mL of PBS media, in a 1.5 mL Eppendorf tube. Due to the lack of a surfactant, the particles cluster together forming ~1-2 µm sized clusters with a concentration of 50 mg/mL. The solution was pipetted thoroughly and subsequently sonicated for 5 min to break apart any aggregation. 2 µL of the suspended clusters is put into 1 mL of E8-Flex medium. Next fluorescent particles (FluoSpheres Carboxylate-Modified Microspheres, 200 nm, Invitrogen) were added to a final concentration of 1:5000 v/v and the solution sonicated for 2 min. The FP tagged the iron oxide clusters forming fluorescent clusters (FCs).

To label the surface of monolayer, immediately after sonication, 100 µL of the FCs were added to each well of a 96 well-plate and incubated for 5 min at 37°C and 5% CO₂. Next the medium was discarded and replaced with E8-Flex medium supplemented with Matrigel at the desired concentration. Control conditions received Matrigel-free medium.

**Microscopy.** A Nikon Ti2 Eclipse inverted microscope with incubation capability was used to image FCs dynamics in various Matrigel conditions and in control conditions using a x10 air objective with 15 min imaging interval for 24 h. Confocal representative images were obtained using a Leica SP8 DIVE (Leica Microsystems) confocal microscope using a x10 air objective or a x20 water-emersion objective. The same microscope was with incubation and resonance scanner capability to capture FP and laminin dynamics using a x10 air objective with 15 min imaging interval for 24 or 48 h. To capture actin
processes and dynamics the same microscope was used with x64 oil-emersion objective with 20 s imaging interval. Low magnification imaging was done using an EVOS (Invitrogen) inverted microscope using a x4 air objective.

**Scanning electron microscopy (SEM).** The media compositions and exposure times followed the same protocol as the pluripotent conditions. First, hPSCs were grown on Matrigel-coated 13mm coverslips inside wells of a 24 well-plate for 24 h under ROCKi conditions, then an additional 24 h under ROCKi-free conditions. At the Matrigel addition timepoint, the Matrigel conditions received E8-Flex medium supplemented with 10% Matrigel and the controls remained Matrigel-free. All media volumes were adjusted to 500 µL due to the well size of the 24 well-plate. After 24 h of Matrigel exposure, cells were fixed with 2.5% glutaraldehyde in 0.1 M Na-cacodylate buffer Ph 7.2. After 3x washing in the same buffer, the coverslips with cells were incubated in 2% osmium tetroxide (Electron Microscopy Sciences, Hatfield PA, US) for 2 h, washed in milli-Q water and dehydrated in a graded ethanol series to 100% ethanol. Finally, keeping the coverslips submerged in 100% ethanol at all times, they were inserted in a coverslip-holder for critical point drying in a Leica CPD300 apparatus (Leica, Vienna, AT). Finally, the dried coverslips were mounted on 25mm SEM support stubs with Carbon-stickers, the edges were silver painted, and the assembly was coated with 10 nm Chromium in a Leica ACE600 coating machine (Leica, Vienna, AT). Cells were studied and imaged in a Zeiss Sigma SEM (Zeiss, Oberkochen, DE) at an accelerating voltage of 2kV and a working distance of 7mm.

**Particle image velocimetry.** PIVlab³⁷,³⁸, was used in combination with MATLAB (MATLAB, R2018b, The MathWorks Inc.). Here FP time lapse images at the Matrigel-monolayer interface were uploaded to PIVlab through graphical user interface. Standard parameter settings were used along with the PIV FFT window deformation algorithm. A 3 pass analysis was implemented and a 128, 64 and 32 pixels interrogation areas were used for Pass 1, 2 and 3 respectively. Calibration was set and later verified using image scale to ensure correct velocity reporting. Data visualization used data smoothing at
maximum strength and the standard strength for the high-pass vector field filter. This was visually
verified to ensure vectors did not misreport velocity orientation and magnitude. Text data extracting
instantaneous and average velocities were obtained. These text matrices were input into MATLAB to
plot the streamlines for the average velocities and vectors for the instantaneous velocities on a contour
map.

**Image analysis.** ImageJ (Version 1.52a) was used for all image analysis. To obtain the velocity dynamics, gray scale velocity maps at 0, 6, 12, 18 and 24 h were extracted from PIVlab and input into ImageJ (Version 1.52a). Radial intensity profiles were obtained using the radial profile angle plugin around domes. The intensities were evaluated for the various time points. 5 velocity zones were radially marked and the velocities within were averaged to obtain 5 velocity values for each time point.

To visualize the trajectories of FCs, 100 FCs per condition were manually tracked using the
MTrackJ plugin on ImageJ. Using the data extracted from the trajectory analysis, average trajectory
velocities and start to finish distances were calculated for each trajectory in different conditions.

Using the OrientationJ<sup>39,40</sup> plugging, ZO1 orientation and coherence scores were obtained.
Using the radial profile angle plugin around domes, radial profiles of ZO1 intensity, FC intensity and
coherence scores were measured for each image in the sequence of images comprising a 24 h time-
lapse. Values from different domes were combined, averaged and compiled into kymographs for ZO1
intensity, FC intensity and coherence scores.

Inter-dome distance evaluation was done by marking the center of individual domes. The
coordinates of each dome were then used to compute the distance to every other dome within the
same field of view. For each dome only the shortest distance is reported, i.e. the nearest dome.
Duplicate values were excluded from the data. A similar technique was used for *in silico* dome
measurement, however here *in silico* coordinates where compared to the *in vitro* coordinates and
nearest dome x-coordinate rank was reported against the *in silico* x-coordinate rank. The average
difference in distance between *in silico* and *in vitro* nearest domes was then evaluated and reported.
To evaluate all radial profiles, using ImageJ, the radial profile function with a sweep angle of 20 degrees and a radius of 100 µm was used with the dome center as the center of the radial profile. The profiles were collected, averaged and reported as a function of radial distance away from the dome center.

To evaluate NANOG or T expressions, using ImageJ, the Hoechst signal was used to segment the nuclei. Each nuclei was rendered an independent region of interest (ROI) and the intensity of the transcription factor expression was obtained and normalized to the Hoechst signal and reported on a scatter plot.

To evaluate the integrated signal of T, confocal z-stacks were summed, the nuclei segmented and the summed intensities of T was reported to compare between Matrigel and control conditions.

**Reaction-diffusion simulation.** A reaction-diffusion (RD) system was built based on the following Brusselator equations for the signals E and P,

$$ \frac{\partial E}{\partial t} = D_E \nabla^2 E + a - E(b+1) + PE^2(1) $$

$$ \frac{\partial P}{\partial t} = D_P \nabla^2 P + Eb - PE^2 (2) $$

The parameters a and b control patterning outcomes from domes to stripes. We chose values of 5 and 7.75 for a and b respectively, which produced periodic domes in continuous domain spaces, and therefore reflecting *in vitro* observations. The E and P diffusivities were set for $D_E = 20$ and $D_P = 160$ respectively which reflected measures *in vitro* inter-dome distances under normal conditions. The diffusivity $D_E$ was changed accordingly to reflected blebbistatin treatment. MATLAB was used to solve a set of Laplacian operations in an iterative process. A spatial increment of 1 unit (10 µm) in both x and y-axis, a timestep of 0.001 s and a total time of 300 s were chosen as they provided a stable simulation. The domain space is by default continuous but can be rendered bounded in any direction by forcing desired edge values to 0 in every iteration. The domain space can be defined as matrix of m rows and n columns but can also be of any arbitrary shape by providing MATLAB a text matrix of 0s, unused spaces, and 1s used domain spaces. This binary matrix can also be obtained through a tissue mask
created in ImageJ using the a saturated Hoescht signal from *in vitro* data. In this way, any tissue shape
can be simulated in the RD model.

**Immunohistochemistry, actin labeling and patterning studies.** After 24 or 48 h cells (pluripotent or mesoderm differentiation conditions) were fixed using 4% paraformaldehyde (Sigma-Aldrich) for 2 h followed by washing with PBS. Sample permeabilization and blocking was made using solution made with 0.3% Triton X (PanREAC AppliChem) and 0.5% BSA (Sigma-Aldrich) for 30 min. Primary antibodies (EZRIN Abcam ab133297 dilution 1:200, NANOG R&D Systems AF1997 dilution 1:200, Brachyury (T) Bioscience X1AO2 dilution 1:200, SOX2 Sigma-Aldrich AB5603 dilution 1:200, SMAD2/3 Bioscience 610842 dilution 1:200, Phosphorylated (p)SMAD1/5/9 Cell Signaling Technologies 13820S dilution 1:200, Laminin Thermofisher PA1-16730 dilution 1:200, 67kDa laminin receptor Abcam ab268079 dilution 1:200) were suspended in permeabilization and blocking solution and applied to the sample for 24 h to ensure diffusion through Matrigel. This was followed by three PBS washes over 24 h. Next immunolabeling was performed over a duration of 24 h using donkey-anti mouse Alexa Fluor 555/647 (Invitrogen dilution 1:500), donkey-anti rabbit Alexa Fluor 555/647 (Invitrogen dilution 1:500), or donkey-anti goat Alexa Fluor 555/647 (Invitrogen dilution 1:500) secondary antibodies. To visualize filamentous actin Alexa Fluor 647 conjugated phalloidin was used (Invitrogen dilution 1:4000). Hoechst was used to visualize nuclei. The ZO1 reported cell line was used to visualize ZO1 expression.

Actin was labeled to visualize live actin dynamics. At 24 h after Matrigel introduction, 100 µL of E8-Flex supplemented with 1 nM of SiR-actin (Spirochrom) was added to the cells and incubated for 10 min. Confocal microscopy was then conducted to observe apical actin dynamics.

The expression of NANOG was used to evaluate regions of high pluripotency in Matrigel and control cases. NANOG and T were used to evaluate and contrast mesoderm expression in differentiation experiments.
**Single-cell RNA sequencing sample preparation.** Samples were grown following the protocol for control and 10% Matrigel conditions in 6 well plates, 3 wells per condition. After 24 h of Matrigel exposure, samples were washed with PBS to ensure little to no Matrigel is left, then dissociated into single cells using an incubation period of 3 min in 1 mL TrypLE Express in every well. The 1 mL solutions were neutralized with 9 mL of DMEM medium supplemented with 20% FBS, the control and Matrigel conditions were pooled, then centrifuged for 5 min at 300 RCF. After decanting the supernatant, the cell pellets were resuspended in 200 µL of E8-Flex medium and put on ice. Library preparations for the scRNAseq was performed using 10X Genomics Chromium Single Cell 3' Kit, v3 (10X Genomics, Pleasanton, CA, USA). The cell count and the viability of the samples were assessed using LUNA dual fluorescence cell counter (Logos Biosystems) and a targeted cell recovery of 10,000 cells was aimed for all the samples. Post cell count and QC, the samples were immediately loaded onto the HyDrop Hybrid⁴¹ platform. Single-cell RNA-seq libraries were prepared using manufacturers recommendations (Single-cell 3' reagent kits v3 user guide; CG00052 Rev B), and at the different check points, the library quality was assessed using Qubit (ThermoFisher) and Bioanalyzer (Agilent). Sequencing was performed using BGI Genomics. The data generated from sequencing was demultiplexed and mapped against human genome reference using CellRanger v3.0.2.

**Single-cell sequencing data processing and clustering:** After QC steps we have retained 6503 and 6746 sequenced cells for the control and 10% Matrigel conditions respectively. We implemented the Seurat⁴² 4.2.0 pipeline to perform scRNAseq analysis. Data normalization was performed, followed by the identification of 2000 highly variable genes using the FindVariableFeatures. Auto scaling of the data was performed and described using principal component analysis (PCA) using the RunPCA function. Data clustering was done using the top 15 PCs and the FindNeighbors function parameter k.param = 20. The FindCluster resolution was set to 0.5. Data visualization was done using the Uniform Manifold Approximation and Projection (UMAP) dimensionality reduction technique using the RunUMAP package with the top 15 PCs. Cluster annotation was performed primarily based on the top 25 genes
and/or specific marker genes, which identified seven clusters, 1) Proliferation, glycolysis, 2) Non-proliferative, Glycolysis, 3) Non-proliferative, 4) Cycling, 5) Transition 1, 6) Transition 2, 7) Neural. Genesets were created for Epiblast, Primitive streak, and Nascent mesoderm identities following the top 25 markers of each of these clusters in the human gastrula dataset²¹.

The 10% Matrigel condition was subset from the scRNAseq dataset, and the data was reanalyzed to obtain a new UMAP using the same parameters as with the full dataset. High level clustering was achieved based on cell cycle markers and NANOG expression, which was shown to be upregulated in domes. Three highlevel clusters were identified expressing high, medium and low nano expressions. Only high level NANOG was considered to be a dome identity, while the others 2 clusters were considered as monolayers.

**Differential gene expression analysis (DGEA), gene ontology (GO) enrichment analysis.** We perform DGEA using the FindMarkers function in the Seurat R packages between 10% Matrigel and control cells with function parameters min.pct = 0.0 and logfc(threshold = 0.0.

The resultant total gene list was then ranked in decreasing expression levels. GO enrichment analyses was then conducted using the GOrilla<sup>43</sup> web portal using the Homo sapiens organism single ranked file running mode and p-value threshold of 10<sup>-5</sup>. The results were then passed to the REVIGO<sup>44</sup> web portal to obtain a summarized GO term list. Default settings were used.

**Single R analysis.** The SingleR<sup>45</sup> package was used to project the cell of the dataset generated here onto the annotated clusters of the human gastrula (E-MTAB-9388)<sup>21</sup>. The top-50 markers genes (log fold change ≥0.5, min.pct ≥0.5) of the clusters identified in human gastrula dataset were used as reference.

**Pseudotime trajectories.** Using the 10% Matrigel subset dataset, the resultant Seurat object was passed to the Slingshot<sup>46</sup> package on R, which was used to obtain and analyze trajectories. Using the
Seurat UMAP coordinates, the Slingshot kmeans function with the parameter centers = 5, the getLineages and getCurves functions, a pseudotime value was obtained for each cell. These values were added to the 10% Matrigel metadata file.

**Quantification and statistical analysis.** Two-way ANOVA statistical tests and unpaired two-tailed t-test with corrections were used where appropriate with a 95% confidence interval (GraphPad Prism 6, Version 6.01, GraphPad Software, Inc.). Statistical significance was considered for all comparisons with *p* < 0.05.

**Data available.** All raw sequencing data, and the combined processed and metadata files generated in this study will become available at GEO.
References

1. Walma, D. A. C. & Yamada, K. M. The extracellular matrix in development. *Development* **147**, (2020).

2. Loganathan, R. et al. Extracellular matrix motion and early morphogenesis. *Dev.* **143**, 2056–2065 (2016).

3. Bénazéraf, B. et al. A random cell motility gradient downstream of FGF controls elongation of an amniote embryo. *Nature* **466**, 248–252 (2010).

4. Palmquist, K. H. et al. Reciprocal cell-ECM dynamics generate supracellular fluidity underlying spontaneous follicle patterning. *Cell* **185**, 1960-1973.e11 (2022).

5. Heisenberg, C. P. & Bellaïche, Y. Forces in tissue morphogenesis and patterning. *Cell* **153**, 948 (2013).

6. Zamir, E. A., Rongish, B. J. & Little, C. D. The ECM moves during primitive streak formation - Computation of ECM versus cellular motion. *PLoS Biol.* **6**, 2163–2171 (2008).

7. Harrisson, F., Van Hoof, J., Vanroelen, C. & Vakaet, L. Transfer of extracellular matrix components between germ layers in chimaeric chicken-quail blastoderms. *Cell Tissue Res.* **239**, 643–649 (1985).

8. Abdel Fattah, A. R. et al. Actuation enhances patterning in human neural tube organoids. *Nat. Commun.* **12**, 1–13 (2021).

9. Gjorevski, N. et al. Designer matrices for intestinal stem cell and organoid culture. *Nature* **539**, 560–564 (2016).

10. Zheng, Y. et al. Dorsal-ventral patterned neural cyst from human pluripotent stem cells in a neurogenic niche. *Sci. Adv.* **5**, 1–14 (2019).

11. Veenvliet, J. V. et al. Mouse embryonic stem cells self-organize into trunk-like structures with neural tube and somites. *Science* (80-.). **370**, (2020).

12. van den Brink, S. C. et al. Single-cell and spatial transcriptomics reveal somitogenesis in gastruloids. *Nature* **582**, (2020).

13. Sanaki-Matsumiya, M. et al. Periodic formation of epithelial somites from human pluripotent stem cells. *Nat. Commun.* **13**, 2325 (2022).

14. Karzbrun, E. et al. Human neural tube morphogenesis in vitro by geometric constraints. *Nature* **599**, 268–272 (2021).

15. Miura, S., Sato, K., Kato-Negishi, M., Teshima, T. & Takeuchi, S. Fluid shear triggers microvilli formation via mechanosensitive activation of TRPV6. *Nat. Commun.* **6**, (2015).

16. Crawley, S. W., Mooseker, M. S. & Tyska, M. J. Shaping the intestinal brush border. *J. Cell Biol.* **207**, 441–451 (2014).

17. Pelaseyed, T. & Bretscher, A. Regulation of actin-based apical structures on epithelial cells. *J. Cell Sci.* **131**, 1–10 (2018).

18. Klingner, C. et al. Isotropic actomyosin dynamics promote organization of the apical cell cortex in epithelial cells. *J. Cell Biol.* **207**, 107–121 (2014).

19. Meenderink, L. M. et al. Actin Dynamics Drive Microvillar Motility and Clustering during Brush Border Assembly. *Dev. Cell* **50**, 545-556.e4 (2019).

20. Akhtar, N. & Streuli, C. H. An integrin-ILK-microtubule network orients cell polarity and lumen formation in glandular epithelium. *Nat. Cell Biol.* **15**, 17–27 (2013).

21. Tyser, R. C. V. et al. Single-cell transcriptomic characterization of a gastrulating human embryo. *Nature* **600**, 285–289 (2021).

22. Soh, G. H., Pomreinke, A. P. & Müller, P. Integration of Nodal and BMP Signaling by Mutual Signaling Effector Antagonism. *Cell Rep.* **31**, (2020).

23. Chhabra, S., Liu, L., Goh, R., Kong, X. & Warmflash, A. Dissecting the dynamics of signaling events in the BMP, WNT, and NODAL cascade during self-organized fate patterning in human gastruloids. *PLoS Biology* **17**, (2019).

24. Streit, A. et al. Chordin regulates primitive streak development and the stability of induced neural cells, but is not sufficient for neural induction in the chick embryo. *Development* **125**, 507–519 (1998).
25. Conlon, F. L. et al. A primary requirement for nodal in the formation and maintenance of the primitive streak in the mouse. **1928**, 1919–1928 (1994).

26. Takenaga, M., Fukumoto, M. & Hori, Y. Regulated nodal signaling promotes differentiation of the definitive endoderm and mesoderm from ES cells. *J. Cell Sci.* **120**, 2078–2090 (2007).

27. Heemskerk, I. et al. Rapid changes in morphogen concentration control self-organized patterning in human embryonic stem cells. *Elife* **8**, 1–28 (2019).

28. Blassberg, R. et al. Sox2 levels regulate the chromatin occupancy of WNT mediators in epiblast progenitors responsible for vertebrate body formation. *Nat. Cell Biol.* **24**, 633–644 (2022).

29. Hart, A. H., Hartley, L., Ibrahim, M. & Robb, L. Identification, Cloning and Expression Analysis of the Pluripotency Promoting Nanog Genes in Mouse and Human. *Dev. Dyn.* **230**, 187–198 (2004).

30. Morkel, M. et al. β-Catenin regulates Cripto- and Wnt3-dependent gene expression programs in mouse axis and mesoderm formation. *Development* **130**, 6283–6294 (2003).

31. Bally-Cuif, L., Gulisano, M., Broccoli, V. & Boncinelli, E. C-Otx2 Is Expressed in Two Different Phases of Gastrulation and Is Sensitive To Retinoic Acid Treatment in Chick Embryo. *Mech. Dev.* **49**, 49–63 (1995).

32. Stavish, D. et al. Generation and trapping of a mesoderm biased state of human pluripotency. *Nat. Commun.* **1**, 1–16 (2020).

33. Bernardo, A. S. et al. BRACHYURY and CDX2 mediate BMP-induced differentiation of human and mouse pluripotent stem cells into embryonic and extraembryonic lineages. *Cell Stem Cell* **9**, 144–155 (2011).

34. Zheng, Y. et al. Controlled modelling of human epiblast and amnion development using stem cells. *Nature* **573**, 421–425 (2019).

35. Muncie, J. M. et al. Mechanical Tension Promotes Formation of Gastrulation-like Nodes and Patterns Mesoderm Specification in Human Embryonic Stem Cells. *Dev. Cell* **55**, 679-694.e11 (2020).

36. Araya, C., Carmona-Fontaine, C. & Clarke, J. D. W. Extracellular matrix couples the convergence movements of mesoderm and neural plate during the early stages of neurulation. *Dev. Dyn.* **245**, 580–589 (2016).

37. Thielicke, W. & Stamhuis, E. J. PIVlab – Towards User-friendly, Affordable and Accurate Digital Particle Image Velocimetry in MATLAB. *J. Open Res. Softw.* **2**, (2014).

38. Thielicke, W. & Sonntag, R. Particle Image Velocimetry for MATLAB: Accuracy and enhanced algorithms in PIVlab. *J. Open Res. Softw.* **9**, 1–14 (2021).

39. Püspöki, Z., Storath, M., Sage, D. & Unser, M. Transforms and operators for directional bioimage analysis: A survey. *Adv. Anat. Embryol. Cell Biol.* **219**, 69–93 (2016).

40. Rezakhaniha, R. et al. Experimental investigation of collagen waviness and orientation in the arterial adventitia using confocal laser scanning microscopy. *Biomech. Model. Mechanobiol.* **11**, 461–473 (2012).

41. De Rop, F. V. et al. HyDrop enables droplet based single-cell ATAC-seq and single-cell RNA-seq using dissolvable hydrogel beads. *Elife* **11**, 1–26 (2022).

42. Satija, R., Farrell, J. A., Gennert, D., Schier, A. F. & Regev, A. Spatial reconstruction of single-cell gene expression data. **33**, 495–502 (2015).

43. Eden, E., Navon, R., Steinfeld, I., Lipson, D. & Yakhini, Z. GOrilla: a tool for discovery and visualization of enriched GO terms in ranked gene lists. *BMC Bioinformatics* **10**, 48 (2009).

44. Supek, F., Bošnjak, M., Škunca, N. & Šmuc, T. REVIGO Summarizes and Visualizes Long Lists of Gene Ontology Terms. *PLoS One* **6**, e21800 (2011).

45. Aran, D. et al. Reference-based analysis of lung single-cell sequencing reveals a transitional profibrotic macrophage. *Nat. Immunol.* **20**, 163–172 (2019).

46. Street, K. et al. Slingshot: cell lineage and pseudotime inference for single-cell transcriptomics. *BMC Genomics* **19**, 477 (2018).
Supplementary Files

This is a list of supplementary files associated with this preprint. Click to download.

* Supp movie1 mov

• Supp movie2 mov

• Supp movie3 mov

• Supp movie4 movie

• Supp movie5 mov

* Supp movie6 mov