Self-assembly of tessellated tissue sheets by growth and collision

Tissues do not exist in isolation–they interact with other tissues within and across organs. While cell-cell interactions have been intensely investigated, less is known about tissue-tissue interactions. Here, we studied collisions between monolayer tissues with different geometries, cell densities, and cell types. First, we determine rules for tissue shape changes during binary collisions and describe complex cell migration at tri-tissue boundaries. Next, we demonstrate that genetically identical tissues displace each other based solely on cell density gradients, and present a physical model of tissue interactions that allows us to estimate the bulk modulus of the tissues from collision dynamics. Finally, we introduce TissEllate, a design tool for self-assembling complex tessellations from arrays of many tissues, and we use cell sheet engineering techniques to transfer these composite tissues like cellular films. Overall, our work provides insight into the mechanics of tissue collisions, harnessing them to engineer tissue composites as designable living materials.


Introduction
A biological tissue is a cellular community or, as Virchow wrote in the 19th century(1), "a cell state in which every cell is a citizen". This concept is increasingly apropos as interdisciplinary research pushes deep into the coordinated cell behaviors underlying even 'simple' tissues. Indeed, cell-cell interactions give rise to behaviors such as contact inhibition (2)(3)(4), collective cell migration (5,6), and cell-cycle regulation (7)(8)(9)(10), which underlie physiological functions such as tissue development and healing (11,12), organ size control (13,14), morphogenetic patterning (15), and even pathological processes such as tumor invasion (16,17). In places where tissues meet, the resulting tissue is a living composite material whose properties depend on its constituent tissues. In particular, tissue-tissue interfaces underlie both biological processes such as organ separation and compartmentalization (18,19), as well as biomedical applications such as tissue-mimetic materials (20)(21)(22) and engineered tissue constructs (23)(24)(25). Thus, recent research has focused on the formation and dynamics of tissue-tissue boundaries. For instance, the interplay between repulsive Eph/ephrin and adhesive cadherin cell-cell interactions regulate tissue boundary roughness, stability, and cell fate (26)(27)(28)(29)(30). Furthermore, colliding monolayers with differences in Ras gene expression were able to displace one another (31,32), while epithelial tissue boundaries were found to induce waves of cell deformation and traction long after the tissues had collided (33). Our goal was to harness these fundamental concepts to define broad 'design principles' for assembling composite tissues in a controlled way. Specifically, we sought to harness mechanical tissue behaviors in the context of cell-sheet engineering, which aims at harvesting intact cell monolayers to create scaffold-free, high-density tissues (24). Such cell sheets are typically produced by allowing cells to come to confluence within a stencil or patterned substrate to form a monolayer with a desired geometry (34,35). Here, we propose an alternative approach where we create arrays of individual epithelial monolayers and then allow them to grow out and collide, fuse at the interfaces, and ultimately selfassemble into tessellated patterns. We performed live imaging as these tissue arrays selfassembled into patterns over 2-3 days, which we predicted by extending our earlier model of tissue expansion (10) to account for multi-tissue interactions. We then characterized the dynamics of the boundary in collisions of tissues with different size, cell density, and composition. Moreover, we proposed a physical model for understanding the resulting boundary motion and extracting tissue mechanical properties from it. We next introduced a design framework for the systematic assembly of many-tissue composites (3 cell types and 30+ tissues), and finally investigated more complex cases such as tissue engulfment and the singular dynamics of tritissue junctions.

Results
Collisions between archetypal tissue pairs. We first characterized interactions between growing pairs of millimeter-scale epithelial tissues, including equal-size rectangles, circles, small vs. large circles, and circles vs. rectangles (Fig. 1, Supplementary Video 1). We used MDCK cells, a standard model (5,7), and labeled each tissue in a pair with a different color (Methods) to clearly distinguish the boundary. Imaging over 2-3 days, we observed no mixing (Supplementary Video 2). Collisions between identical rectangular tissues are well characterized from the traditional wound healing scratch and barrier removal assays (5,33), and our data here confirmed the expected symmetric collision and fusion along the mid- line (Fig. 1a). We next compared identical circles, where we observed a straight boundary form at the midline as before (Fig. 1b), but it was almost two times smoother than the boundary between colliding rectangles ( Supplementary Fig.  1a, p value 0.014, see Methods). We suspect that this is due to the difference in collision dynamics: while parallel strips of tissue collide all along the collision line at once, circles collide at a single point and gradually extend the boundary line outward (Supplementary Videos 2-3). We introduced asymmetry by replacing one of these circles with either a much larger circle, or a long, thin rectangle (Fig. 1c,d). In each case, we observed a curving of the boundary away from the larger tissue, which was especially notable in the circle-rectangle collisions. We aligned and averaged the final segmented fluorescence signals to demonstrate how stereotyped these collision patterns are (Fig. 1, rightmost column).
Predicting the shape of colliding tissues. The stereotyped nature of these collision patterns suggested value in a computational design tool to predict the evolution of tissue shapes upon collisions. We previously established that freely-growing epithelia expand outward with a normal velocity v n , which, except in high-curvature zones, is uniform around the perimeter of a tissue and independent of the tissue geometry or density (10). Here, we implemented this observation into our model to predict the expansion and in-teraction of multiple tissues by assuming that tissue edges pin in place upon contact (Supplementary Fig. 2 and Supplementary Note 1). We initialized the model simulations using the initial tissue locations from experiments, and we used v n = 29.5 µm/h as measured in Ref. (10) and confirmed here (Methods). Without any fit parameters, these simulations predict the shape evolution of the colliding tissue pairs in our experiments ( Fig. 1, Supplementary Video 4, blue/yellow/white outlines show model predictions). Consistent with our observations, pairs of equally-sized rectangles or circles produce a straight boundary, while mismatched shapes produce a curved boundary (Fig. 1a-d). In our model, this is because the initial tissue edges are equidistant to the dividing midline in equal tissues but closer to the midline in large tissues than smaller tissues. In all cases, we found that the mean error of the predicted boundary was compatible with its measured roughness ( Supplementary Fig.  1b), showing that our modeling approach is appropriate at these large scales.

Homotypic tissue boundary dynamics and collision memory.
Having analyzed the macro-scale patterns formed by colliding tissues, we next focused on the dynamics at the collision boundaries. Using the same configuration of two rectangles as in Fig. 1a as a control (Fig. 2a), we then compared the boundary dynamics of tissue pairs with a mismatch in either tissue width (500 µm vs. 1000 µm, Fig. 2b) or cell density (2640 cells/mm 2 vs. 1840 cells/mm 2 , Fig. 2c). First, we determined how asymmetry in tissue width or density affected boundary motion upon collision. We tracked the mean tissue boundary and found that wider and denser tissues displaced narrower and less dense tissues, respectively. Boundary motion was pronounced, directed, and sustained for 15−20 h before stopping (red and blue in Fig. 2d; Supplementary Videos 5-6). In contrast, control experiments with symmetric tissue collisions showed larger boundary fluctuations with very little average drift (black in Fig. 2d; Supplementary Video 7). Prior studies have noted similarly biased boundary dynamics, but only in heterotypic tissue collisions, for example between wild-type and Ras-transformed endothelial cells (31,32). Here, we show that collisions between homotypic tissues -genetically identical -also produce boundary motion due to asymmetry in tissue size or cell density. In contrast to heterotypic collisions (32), however, homotypic tissue boundary motion eventually stops. We related boundary motion to tissue flow using particle image velocimetry (PIV) to measure the velocity field. We represented these data in kymographs of the velocity component along the collision direction, v x , averaged over the tissue length, across multiple tissue collisions (Fig. 2e, see Methods). With identical (control) tissues, cells around the tissue boundary symmetrically reversed their velocity shortly after collision; convergent motion became divergent. We defined the "center of expansion" as the position from which tissue flow diverges (Methods). For control tissues, the center of expansion lies very near the tissue centroid shortly after collision (Fig. 2e left).  At what point, if any, do two fused tissues act as one? We investigate this question in collisions between tissues with size or density mismatch, which exhibited tissue flow towards the less dense or narrower tissues. In these cases, the centers of expansion began at the centroid of wider or denser tissues rather than at the overall centroid or collision boundary ( Fig. 2e center and right). The center of expansion then gradually shifted towards the centroid of the fused tissue. After the center of expansion reached the overall centroid, the fused tissue expanded symmetrically without memory of the collision. Thus, by comparing expansion centers to geometric centroids, we identified the transition whereby two colliding tissues shift behaviors to act as one larger tissue.
Cell density gradients drive boundary motion. We hypothesized that collision boundary motion was driven by cell density gradients (36)(37)(38)(39)(40). To test this, we quantified local cell density (Methods) and represented it in kymographs for each collision assay (Fig. 2f). In all cases, we found collision boundaries moved down local density gradients, consistent with our hypothesis. While tissues in the size-mismatch assay had the same initial density, the larger tissue had a higher density at the time of collision ( Fig. 2f center). This observation is consistent with our prior work showing that, even when prepared with the same density, larger tissues develop higher cell densities than smaller tissues as they expand (10).
To understand the mechanics of this process, we modeled the expanding tissue as an active compressible medium (41). Tissue expansion is driven by polarized active cell-substrate forces, which are known to be maximal at the tissue edge and decay over a distance L c ∼ 50 µm as we move into the cell monolayer (42,43). Hence, we ignore active traction forces at the tissue boundary after collision, which is ∼ 1 mm away from the outer tissue edges. At the collision boundary, we establish a force balance whereby pressure gradients drive tissue flow v as where ξ is the cell-substrate friction coefficient. Moreover, we assume that the tissue pressure P increases with cell density ρ as specified by an unknown equation of state P (ρ), which predicts that the collision boundary moves from high to low densities with a speed proportional to the density gradient (Fig. 2g).
To test this prediction, we measure both the velocity and the density gradient at the boundary for each experiment in our three different assays (Methods). Consistent with our prediction, the results show a negative correlation between the boundary velocity and the cell density gradient (Fig. 2h).
Estimating tissue mechanical properties from collisions. Based on our model, we use our measurements of cell density and velocity to extract information about the tissue's equation of state P (ρ). To this end, we obtain the average boundary velocity and density gradient for each assay, and we fit a line to them (Fig. 2h). From this fit, and using ξ ∼ 100 Pa·s/µm 2 (43, 44), we obtain P (ρ) ∼ 1.5 Pa·mm 2 . This result indicates that, in the conditions of our experiments, for every cell that we add per square millimeter, the tissue pressure goes up about 1.5 Pascals. Next, we use these results to estimate the mechanical properties of the cell monolayer. To this end, we assume a specific equation of state: where K is the bulk modulus of the monolayer and ρ e is a reference cell density. This equation of state was justified theoretically for growing tissues around their homeostatic state, around which the cell proliferation rate varies linearly with cell density (45    compression that results from a tissue collision. Overall, our collision assays provide a way to probe the bulk mechanical properties of migrating cell monolayers, which are otherwise difficult to measure. Remarkably, analyzing collisions between tissues that differ only in their cell density allowed us to infer mechanical properties without measuring any mechanical forces. Rather, we employed our model to relate tissue flows to pressure and density gradients, from which we inferred the relationship between pressure and density. In the future, collision assays might be used to measure the equation of state of cell monolayers, which is a key input for mechanical models of growing and expanding tissues (41,49).
Heterotypic tissue boundary dynamics. Having studied collisions between homotypic tissues, we now turn to collisions between heterotypic tissues with different cell migration speed and cell-cell adhesion strength. We prepared co-cultures of the breast cancer cell lines MCF10A (benign), MCF7 (malignant, non-invasive), and MDA-MB-231 (metastatic) as monolayers of the same size and cell density. We used homotypic MCF10A collisions as a reference, for which we observed non-mixing and boundary dynamics similar to the homotypic MDCK collisions discussed earlier (Supplementary Video 8).
We first collided monolayers of MCF10A and MDA-MB-231 cells, which have the largest phenotypic difference among the cell lines we used. While these tissues have similar expansion speeds, they exhibit radically different collective dynamics, reflective of different cell-cell adhesion strengths (50) (Supplementary Video 9). While cells in MCF10A tissues hardly exchange neighbors, the metastatic MDA-MB-231 cells con-tinually undergo neighbor exchanges and even crawl over each other out of plane. Upon collision, the MCF10A tissue simultaneously displaced and crawled underneath the MDA-MB-231 tissue (Fig. 3a-d).
We next investigated collisions between MCF10A and MCF7 monolayers. The MCF7 monolayer expands about 6 times more slowly than the MCF10 monolayer, which allowed us to explore the effects of different edge speeds on tissue collisions. Surprisingly, we found that the slower MCF7 tissues actually displaced the MCF10A tissues (Fig. 3e, Supplementary Video 10), which may be due to differences in cell-cell adhesion and eph/Ephrin signaling (51). This, at least, shows that a higher expansion speed does not imply a higher "strength" upon collision. In fact, MCF10A cells at the collision boundary reversed their velocity and migrated away from the MCF7 tissue within 8 hours after collision, starting at the boundary and progressively moving into the MCF10A monolayer (Fig. 3f). This behavior seems a tissuescale analog of the cellular behavior known as contact inhibition of locomotion, whereby a cell stops and reverses its direction of motion upon collision with another cell (2)(3)(4).
Further, in collisions between tissues with different expansion speed, the faster tissue should be able to engulf the slower tissue, similar to the engulfment between tissues with differential adhesion (52). We confirmed this hypothesis in collisions between strips of MCF10A cells (fast) and circles of MCF7 cells (slow), which we reproduced with our tissue shape model by incorporating different speeds into our simulation ( Fig. 3g-  tissue interfaces, but here we highlight how differences in expansion speed enable design options. Large-scale tissue tessellations for cell sheet engineering. The stereotyped nature of tissue-tissue collisions suggest simple underlying design rules that would allow selfassembled tissue tessellations to be designed first in silico and then realized in vitro. We tested this idea with a tesselation inspired by the artwork of M.C. Escher-a 'dice lattice' (Fig. 4a). To design this tesselation, we used the computational model described above to simulate many initial tissue array conditions until converging on the pattern of ellipses shown in Fig. 4b. We engineered this pattern with tissues ( Fig. 4c) and filmed it developing as predicted ( Fig. 4d; Supplementary Video 12). This computer-aided-design (CAD) process generalized to arbitrary tessellations (Fig. 4e,f; Supplementary Videos 13-14), offering a 'TissueCAD' approach to designing and building composite tissues. Composite cell sheets may be particularly useful for tissue engineering, where cell sheets are extracted from culture vessels and used as building blocks for larger constructs. We demonstrated compatibility of this process with our tissue composites by culturing a dice lattice on a temperatureresponsive substrate (UpCell dishes) and then transferring the tissue to a new culture dish (Fig. 5, Methods). The morphology of sharp tissue-tissue interfaces were preserved during the transfer, demonstrating that such tissue composites can, in principle, be handled like standard cell sheets.
Dynamics at tri-tissue collisions. During the selfassembly of tissue tesselations, we observed a special behavior at tri-tissue collision points: We often found long streams of tissue that necked down to the single cell scale, visually reminiscent of streams of invasive cancer cells (Fig. 6a-b). However, these events, which we call 'escapes', involve the co-migration of all three tissues rather than the invasion of one tissue into the others (Fig. 6c-e, Supplementary Videos  15-16). Consequently, the initial relative positions of the three colliding tissues is a strong statistical determinant of escape events (Supplementary Fig. 5). We characterized the dynamics of escapes by measuring cell speed fields around tri-tissue collisions, which showed that the escaping tissue migrated faster than its neighbors ( Fig. 6fh, see Methods). To determine whether this speed increase was a generic consequence of the local negative curvature of tri-tissue collision points, we compared tri-tissue collisions to a single tissue patterned to match the overall shape of the colliding tissues at the time of escape (Fig. 6i). For the single tissue case, we did not observe any speed increase (Fig. 6h,j), which rules out local curvature as the sole cause of escape events. These findings suggest that escapes are an emergent dynamical property of three-tissue interactions.
Overall, multi-tissue collisions produce unique, dynamic boundary conditions and mechanical states that give rise to surprising, almost morphogenetic behaviors at tissue junctions, suggesting an important role for collisions in composite tissue development and engineering.

Discussion
We investigated how tissue-tissue interactions can be harnessed to self-assemble complex composite tissue sheetstissue tessellations. First, we demonstrated that colliding tissues change shape in stereotyped and predictable ways. Then, we proposed a physical model of tissue-tissue collisions that links the motion of the tissue boundary to underlying gradients in cell density, which drive tissue flow. Further, we used this model to estimate material properties of the colliding tissues without any force measurements. In this con-text, previous work had shown that heterotypic tissues can displace each other upon collision (31,32). Our findings revealed that even homotypic tissues, which are genetically identical, can displace each other based on purely mechanical differences. Therefore, our collision assays could be used to study mechanical tissue competition (36-40, 53, 54), which might provide biophysical insight into development (55,56), homeostasis (57), and tumor growth (36,53). Based on the reproducible and almost algorithmic tissue interactions that we found, we developed computational design tools to create complex tissue tesselations. The tesselations are obtained by self-assembly, which allows the tissue boundaries to develop naturally. Thus, our work demonstrates how tissue sheets can be treated as 'designable' living materials. Specifically, we developed a simulator that, despite lacking both biophysical laws and cellular resolution, predicts tissue patterns accurately at the 100+ µm scale. This feature makes the simulator useful to design tissue composites in silico before realizing them in vitro. This approach is compatible with advanced biofabrication strategies such as cell sheet engineering, which we demonstrated by transferring an 'Escher' tissue sheet between Petri dishes while preserving tissue integrity and internal boundaries. Tissue tessellation should also be compatible with bioprinting, which could be used to pattern larger arrays of the initial tissue seeds.

Materials and Methods
Cell culture. MDCK.2 wild type canine kidney epithelial cells (ATCC) were cultured in customized media consisting of low-glucose (1 g/L) DMEM with phenol red (Gibco, USA), 1 g/L sodium bicarbonate, 1% streptomycin/penicillin (Gibco, USA), and 10% fetal bovine serum (Atlanta Biological, USA). MCF10A human mammary epithelial cells (ATCC) were cultured in 1:1 DMEM/F-12 (Thermo Fisher Scientific, USA) media which consists of 2.50 mM L-Glutamine and 15 mM HEPES buffer. This media was supplemented with 5% horse serum (Gibco, New Zealand origin), 20 ng/mL human EGF (Sigma, USA), 0.5 µg/mL hydrocortisone (Fisher Scientific), 100 ng/mL cholera toxin (Sigma), 10 µg/mL insulin (Sigma, USA), and 1% penicillin/streptomycin (Gibco, USA). MDA-MB-231 (ATCC) and MCF7 (ATCC) human mammary cancer cells were both cultured in 1:1 DMEM/F-12 (Thermo Fisher Scientific, USA) media supplemented with 10% fetal bovine serum (Atlanta Biological, USA) and 1% penicillin/streptomycin (Gibco, USA). All cells were maintained at 37°C and 5% CO 2 in humidified air. Tissue dye segmentation. The tissue dye becomes diluted as cells divide and spread, so the dye at tissue edges (where cells are more spread and divide more frequently) becomes much more dim than the center of tissues. Because we saw no mixing in our collisions, we segmented the fluorescence channels using a custom MATLAB (Mathworks) script and overlaid them with the phase contrast images for clear visualization. To segment fluorescence images, we normalized the fluorescence channel histograms to each other and compared relative brightness for each pixel between channels. We then masked with the binary masks obtained from the phase contrast channel.
Setting v n for model. We set the normal velocity for the model (all shapes and tessellations) according to the outward velocity of outer edges of the control rectangle collisions. The outward velocity was found to be 29.4 ± 2.3 µm/h (standard deviation), so we used the previously reported speed for expanding circles of 29.5 µm/h (10).
Velocity measurements. We calculated tissue velocity vector fields from phase contrast image sequences, rotating each image so that the initial tissue locations in image pairs were horizontal. We used the free MATLAB package PIVLab with the FFT window deformation algorithm, employing a 1st pass window size of 96×96 pixels and second pass of 48×48 pixels, with 50% pixel overlaps. This resulted in a final window size of 88×88 µm. Data was smoothed in time with a moving average of 3 time points centered at each timepoint.
Average kymographs. We first constructed kymographs of each rectangular collision pair, averaging over the vertical direction of each timepoint and ignoring the top and bottom 1 mm. We then averaged the individual tissue kymographs, aligning by the initial tissue configuration, and determined the edge extent from the median extent of the individual kymographs.

Center of expansion.
We determined the center of expansion by thresholding as |v x | < 3µm/h for individual kymographs of v x . We filtered for the largest contiguous region and took the midline of this region as the center of expansion.
Cell density measurements. We first reproduced nuclei positions from 4X phase contrast images using our in-house Fluorescence Reconstruction Microscopy tool (58). The output of this neural network was then segmented in ImageJ to determine nuclei footprints and centroids. Local density was calculated for each PIV window by counting the number of nucleus centroids in that window.
Boundary velocity and cell density gradient determination. Boundary velocity was found from the position change of the midine in Fig. 2d. Cell density gradient ∂ x ρ was found as ρ R −ρ L x R −x L , where ρ L and ρ R are the total density within 300mum wide regions immediately to the left and right of the tissue boundary, respectively, and x R − x L = 300µm. We plotted ∂ x ρ for timepoints between 20 h and 36 h, which is after collisions and before the boundary stops moving.
Cell sheet engineering tissue patterning and transfer. We first patterned tissues on 3.5 cm NUNC TM UpCell TM dish with supportive membrane (ThermoFisher Scientific, USA). We followed the same collagen coating and stencil application process as before, but passivated the underside of our stencils to avoid damaging the UpCell TM surface. To passivate the stencils, we incubated them for 30 min at 37°C in Pluronic TM F-127 solution (ThermoFisher Scientific, USA) diluted in PBS to 2%. We washed the stencils three times in DI and gently dried them with compressed air before transferring to the dish. After the tissues reached confluence within the stencils, we removed the stencils and allowed the tessellation to collide and heal for ∼72 h. To release the tessellation monolayer from the dish, we changed to cold media and moved the dish to an incubator set to 25°C for 1.5 h. We then pre-soaked the supportive membrane in media to avoiding membrane folding, and floated the membranes on the media above the tessellations. We then carefully aspirated the media from beside the membranes to ensure tight contact between the membrane and monolayer with no bubbles. We moved the UpCell TM dish with membrane to a 4°C refrigerator to ensure total release, and prepared a standard 3.5mm tissue culture dish (BD Falcon, USA) coated with collagen IV as before and filled with warm media. After 7 minutes at 4°C, we then carefully removed the membrane and tessellation monolayer from the UpCell TM dish and floated it in the tissue culture dish with the tessellation side down. We aspirated the media from beside the membrane to initiate bubble-free contact with the dish surface and covered the membrane with 350 µL of warm media. We incubated the membranes at 37°C overnight before floating the membrane off the surface by filling the dish with media and removing it with tweezers.