Collective movement and organization of cell monolayers are important for wound healing and tissue development. Recent experiments highlighted the importance of liquid crystal order within these layers, suggesting that +1 topological defects have a role in organizing tissue morphogenesis. We study fibroblast organization, motion, and proliferation on a substrate with micron-sized ridges that induce +1 and -1 topological defects using simulation and experiment. We model cells as self-propelled deformable ellipses that interact via a Gay-Berne potential. Unlike earlier work on other cell types, we see that density variation near defects is not explained by collective migration. We propose instead that fibroblasts have different division rates depending on their area and aspect ratio. This model captures key features of our previous experiments: the alignment quality worsens at high cell density and, at the center of the +1 defects, cells can adopt either highly anisotropic or primarily isotropic morphologies. Experiments performed with different ridge heights confirm a prediction of this model: Suppressing migration across ridges promotes higher cell density at the +1 defect. Our work enables a mechanism for tissue patterning using topological defects without relying on cell migration.

collective migration | topological defects | cell motility | pattern formation Monolayers of cells in multicellular organisms cooperate to transmit forces in embryogenesis, act as a barrier, and perform many more essential functions (1). These cells often have long axes locally aligned with each other-i.e., they have local nematic order akin to liquid crystals (2)(3)(4)(5). Deviations from perfect nematic alignment can occur as topological defects. In 2D, defects are points where following cell orientation for a complete cycle around the defect leads to a rotation in orientation θ = 2π q, where the topological charge q is integer (q = ±1 shown in Fig. 1A) or half-integer. Topological defects are biologically relevant: They can drive cell death and extrusion (6), cell dynamics (7), tissue branching in regeneration (8), and growth (9). These defects can also reorganize cell density. Recent experiments with monolayers of various cell types show that cells tend to congregate at positive defects and disperse at negative ones (9)(10)(11), though this is not universal to all cell types (6,12). Congregation at positive defects can result in increased density at +1/2 defects (10), creation of new layers of cells at +1/2 defects (11), or growth of mounds of cells at +1 defects (9). In all these examples, accumulation at defects with positive topological charge and depletion at defects with negative topological charge is driven by collective migration of cells. Here, we want to understand how we can control cell density, shape, and cell orientation by exploiting the topological properties of cell monolayers. We use our earlier-developed system of NIH 3T6 fibroblasts on a substrate with micron-scale ridges (12). Fibroblasts are spindle-shaped cells that, apart from steric interactions, rarely interact with each other, not developing strong cell-cell adhesion (4,13). Although highly motile at low densities, fibroblast speeds quickly decay such that at high densities cell arrangement can be characterized by theories of nematic liquid crystals in equilibrium (13). Fibroblasts on 2D substrates exhibit back and forth motion along their long axis without preferential direction (4,13) and are thus significantly less polar than cell types that develop persistent collective migration (9,14). Therefore, we might expect fibroblasts to have qualitatively different responses to induced topological defects than, e.g., myoblasts or epithelial cells.

When we impose a +1 or -1 defect pattern using the ridged substrate (Fig. 1A), fibroblasts align their long axes along the ridges to take up this defect pattern (Fig. 1B). We also see enhanced density at +1 defects and decreased density at -1 defects. However, we find that system lacks large-scale collective migration that could drive density variations at defects. We use simulations and experiments to show that shapedependent division is sufficient to cause density variations at defects. This is a qualitatively different way to pattern cells using topological defects.

When elongated cells exist at high densities, such as in tissues, their long axes tend to align with each other. This order is not perfect, and places where it failsdefects-can be used to control the tissue's properties, guiding cell death and tissue motion. Here, we place elongated fibroblasts on ridged patterns that induce defects. These defects change cell density-but not by changing cell crawling motion, as found in other cell types. We instead argue that fibroblast division is sensitive to cell shape and cell shape is changed by the pattern we use. We show using modeling that this process can explain our experimental results. Our work provides an additional set of tools to control and organize tissues.

Likely Not through Migration. We use photolithography to create a substrate with 1.5-μm-high ridges in a pattern chosen to induce +1 and -1 topological defects (Fig. 1 A and B). We seed 3T6 fibroblast cells on the fibronectin-coated substrates and observe their behavior as they proliferate. Our earlier work (12) discovered three key features, which we reproduce in Fig. 1 C and D (gray box): 1) the fibroblasts' long axes follow the ridges, 2) the degree of deviation from the ridges increases as cells are increasingly packed past confluence, and 3) fibroblast density is increased relative to the rest of the monolayer at the center of the +1 defect pattern and relatively decreased at the center of the -1 defect pattern. It is apparent in Fig. 1B that the long axis of the fibroblasts follows the direction imposed by the ridges-though imperfectly. We measure the quality of cell alignment via the RMSD of cell long axis orientation θ from expected alignment θ e (Materials and Methods), RMSD ≡ δ 2 ≡ (θθ e ) 2 . Cells at subconfluent density ρ < 1,000 cells/mm 2 are relatively wellaligned, but the quality of alignment decreases at larger densities (Fig. 1C).

We also measure cell density, via nuclear fluorescence I (r), as a function of the distance r from the pattern center. Fibroblasts have an elevated density close to the +1 defect and low density at the -1 defect (Fig. 1D).

How does the topographic pattern change density at the defect core? Previous studies on myoblasts, neural progenitor cells, and myxobacteria argued that density differences at defects arise due to collective motion of cells (9)(10)(11), including dramatic inspiraling migration (9). We sketch this broad mechanism in Fig. 1E. However, increased density at the +1 defect could also arise from higher proliferation rate of cells near the +1 defect (Fig. 1F ). If migration were driving the increase in density in our experiments as in other cell types, we would expect significant inward migration toward +1 defects and away from -1 defects. We find instead that 3T6 fibroblasts primarily move azimuthally around the +1 defect, with short-range correlation of velocities, but no broad inward flow (Fig. 1 G and H, SI Appendix, Fig. S1, and Movie S1). Because roughly equal numbers of cells are moving clockwise and counterclockwise and there is no collective inward spiraling motion or relevant net inward motion, we hypothesize that cell division rate differences are the driving factor of our observed density differences, and we develop a model with this assumption.

ment Patterns. We model spindle-shaped fibroblasts as selfpropelled deformable elliptical particles. Cell i has semimajor axis length a i and semiminor axis length b i (Fig. 2A). Our model includes cell motility, cell-cell interactions, and cell division: We give a full description in the Materials and Methods, and a brief summary here. At every step in the model, cells are chosen randomly to either move, rotate, or alter one of their axis lengths (Fig. 2B) and we accept this move with a probability depending on the change in system energy. The energy is determined by cell shape, cell-cell interactions, cell-substrate interactions, and cell polarity. One attempt for each cell is a "Monte Carlo step" (MCS); we calibrate parameters so 100 MCS is 1.5 min of experimental time.

The cell shape energy models cells tending to keep their area A = π ab and aspect ratio AR = a/b; we set AR pref = 4, A pref = 1,400 μm 2 to match experiment (4). Deviations from preferred aspect ratios and areas have energy cost proportional to k AR and k Area , respectively, setting the "stiffness" of the cell's shape. We also include a core energy that prevents indefinite squeezing of cells. Cells interact with one another via a modified Gay-Berne (15-18) potential widely used in liquid crystal simulations (Fig. 2C); this energy promotes cells having their long sides adjacent to one another, inducing nematic order. Cell-ridge overlap is penalized with energy cost equal to the product of ridge strength k r and fraction of the cell overlapping with the ridge. We argue that ridge strength reflects the ridge height in experiments. We have chosen the ridge structure in the simulations to resemble the ridges used to induce +1 and -1 topological defects in experiments.

Crawling eukaryotic cells are chemically and mechanically polarized (19); we summarize this polarity by a vector p indicating the direction the cell prefers to move (Fig. 2A). Fibroblasts move along their long axis (4). We add a motility energy that encourages motion along p and the long axis in a direction = û( û • p), where û is the long axis of the cell. Fibroblast polarity occasionally flips direction ( 4), which we model by stochastically reversing p with average flip period of ∼2.5 h. We also assume cells tend to align to their past displacement (20)(21)(22), leading to some coordination between cell velocities, as in Fig. 1G. In between polarity flips, p obeys a rule proposed by Szabo et al. (23) where after t MCS, we update the polarity vector for each cell as p t = (1 -1/τ pol )p t-1 + r where r = ( x, y) is the proposed displacement, p t-1 is polarity for that cell at time-step t -1 and τ pol is the polarity decay timescale measured in MCS. Polarity thus reorients toward the most recent displacement, promoting cells crawling persistently and coherently (21). This is disrupted by polarity flipping.

Within our model, we seed initially small, circular cells at density of ρ init ∼ 70 cells/mm 2 randomly in our periodic simulation box and then choose cells to divide at a rate set by the experimental growth curve (SI Appendix, Fig. S2). The probability that cell i is selected for division is p i , given by Eq. 1, discussed in more detail in the next section.

We show a typical simulation in Fig. 2 D and E and Movies S2 and S3. The -1 defect is constructed by using the periodic boundary condition (SI Appendix, Fig. S3). We track simulated cells' RMSD from perfect alignment with the ridge pattern δ 2 as cells proliferate and observe that RMSD increases as cells reach higher density (Fig. 3A). This is consistent with our experiments (Fig. 1C). In our model, alignment is also controlled by ridge strength k r -larger k r decreases RMSD (Fig. 3A), though this effect saturates for k r > 100. temperature T ; Materials and Methods.) We get almost identical evolution of RMSDs if we compute cell orientation deviations from the orientation of the closest ridge point (SI Appendix, Fig. S4). The decrease of alignment as cells proliferate and the dependence of alignment on ridge strength are simply a consequence of packing anisotropic deformable cells and do not require cell motility or a particular division mechanism (SI Appendix, Fig. S5). As cells proliferate and alignment worsens, the average aspect ratio of cells decreases (SI Appendix, Fig. S6A) and cells with smaller aspect ratios deviate more from perfect alignment than elongated ones (SI Appendix, Fig. S6C). This is consistent with earlier experimental observations that epithelial cells-which are more isotropic-have higher RMSDs (12). Ridge alignment thus decreases at large densities because cells become less anisotropic-hence less able to coherently align. The experimental analog to ridge strength is ridge height above the substrate. We vary ridge height experimentally (Fig. 3B), finding weak effects on RMSD, suggesting experiments are near the limit where increasing k r has diminishing returns on alignment. Simulated cell motion resembles experimental trajectories (Fig. 1G). We see locally correlated, primarily azimuthal motion without overall coherent direction (Fig. 3 C and D). Average cell speeds slow over time as the monolayer becomes more densely packed, broadly consistent with past measurements (Fig. 3E) (4).

Our model recapitulates experimental cell motion and alignment. Can we understand the increase in density at +1 defects and decrease in density at -1 defects?

Variations at Defects. Higher density of cells near +1 and lower density at -1 defects could arise from cell migration or cell proliferation (Fig. 1 E and F ). Given the lack of clear inward migration (Fig. 1 G and H ), we hypothesize cell proliferation rates are different near defects. One possible reason for this difference is that the cell shapes near the defects differ. Past work on confinement and stretching experiments with endothelial and smooth muscle cells demonstrated that decreased spread area suppresses proliferation while uniaxially extended cells-with large aspect ratio-also have suppressed proliferation (24,25). We thus propose a model where larger cells and more isotropic cells are more likely to divide. We set the probability of cell i to be selected to divide as

where the shape sensitivity α tunes how sensitive division is to cell shape and Z is a normalization factor. When α = 0, cells with the biggest area are the most likely to divide independent of AR; as α increases from 0 to 2, more isotropic cells (AR → 1) with large areas become more likely to divide (Fig. 4A). Since cells are more likely to become isotropic when densely packed within the inner rings of the +1 defect, and isotropic cells more likely to divide, shape-dependent division can potentially drive the density variations. To test whether shape-dependent division is sufficient to reproduce density variation, we vary sensitivity to shape α and observe density changes near defects. We show the change in density relative to the whole-system average density ρ avg in Fig. 4B-analogous to experiments in Fig. 1D. When division probability is independent of shape (α = 0), density does not strongly depend on distance from the +1 defect, but density is below its average value near the -1 defect. (Similar results are found when cells are randomly selected to divide: see SI Appendix, Fig. S7.) When we increase α → 2, making isotropic cells more likely to divide, relative density at the +1 defect center grows significantly. There is also a slight increase in density near the -1 defect for larger α (Fig. 4B), but the normalized density deviation from the average remains negative. Cell shapes also change. As α → 2, we see that cells near the +1 defect become more isotropic than their surrounding cells. On the other hand, cells at the core of the -1 defect are always more elongated than those further away (Fig. 4C). These patterns are consistent with our expectation that, when α = 2, isotropic cells that are more likely to divide are near +1 defects and that decreased aspect ratios allow more cells to pack near +1 defects. This creates a positive feedback loop between cell shape and density, where isotropic packing leads to higher density and higher density leads to more isotropic packing. In general, cells at the core of the +1 defect (which are more isotropic) have higher deviations from their expected alignment while cells near -1 defect (which are more elongated) have better alignment (SI Appendix, Fig. S6B). As cells become more isotropic at higher densities near the +1 defect, the cell-cell interaction potential becomes less dependent on cell orientation. In the limit where cells become circular, they don't have preferred alignment and orient randomly, which we believe causes larger RMSDs. On the other hand, elongated cells near the -1 defect align well with each other and the ridges, resulting in lower deviations from expected orientation.

Increased density at the +1 defect is made more prominent by cell motility, but can also be seen without it (k move = 0; SI Appendix, Fig. S8). In the absence of motility, cells rarely cross ridges, so the ∼40% of simulations that start with no cell in the inner ring still have low density in the inner ring at the simulation end. Even though density at the +1 defect may be high in the other 60% of simulations, the large fraction of simulations with zero or low density at the core means the overall density increase at the +1 defect is weaker in the absence of motility.

Shapes. When α = 2 in our simulation, cells near the defect are more isotropic than those further away (Fig. 4C), but this is highly variable between simulation runs. Increasing α from 0 to 2 increases the number of simulations with isotropic cells at the core-but there are many cases with anisotropic cells (Fig. 5A). We also see similar variability experimentally. In both experiments and simulations, we observe both patterns with elongated cells at the +1 defect center (Fig. 5B) and isotropically packed cells at the +1 defect (Fig. 5C).

Height Enhances Density Variations. How can we test our idea that division is driving density variations? If collective migration is the cause of density variations, we expect preventing cells from crossing ridges to suppress density changes. * On the other hand, if shape-dependent division drives density differences, then density variations should increase if we constrain cell movement across the ridges. At first, in our simulation, we increase ridge strength k r from 60 to 120 to reduce ridge crossing (illustrated in Movies S4 and S5), finding that the number of cells overlapping with ridges drops dramatically as ridge strength grows (Fig. 6 A and B). The decrease in ridge crossing is accompanied with marked changes in density (Fig. 6C), aspect ratio, and alignment (SI Appendix, Fig. S9) near defects. We see a relatively uniform density near the +1 defect for weak ridge strength k r = 60, but we see muchincreased relative density at the +1 defect as k r → 120 (Fig. 6C). We see the opposite trend near -1 defects, with relative density decreasing with k r , though this saturates as k r ≈ 80 -120. Thus, simulations predict that preventing cell crossings enhances density at the core of the +1 defect-as expected if cell division drives the density increase.

To experimentally test this prediction, we vary ridge height to constrain cell movements across ridges. Confocal microscopy of cell nuclei shows a reduction of fibroblast-ridge overlap as we increase ridge height from h = 1.5 μm to h = 14 μm (Fig. 6D). We quantify this by measuring the ratio f on /f off between the * We note that, in principle, increasing ridge height could increase density changes even when driven by motility. This would happen if we were moving from no ridges to a weak ridge height. However, we are not likely to be in a situation where small ridges are a perturbation to no ridges. We see that increasing ridge height does not increase alignment (Fig. 3B), and increasing ridge height decreases crossings (Fig. 6). We thus expect increasing ridge height would suppress density changes if collective motility were driving the density increase at the +1 defect. fraction of the ridge area occupied by nuclei (f on ) and the fraction of nonridge area occupied by nuclei (f off ) (Fig. 6E). The average f on /f off decreases roughly by a factor of four as the ridge height goes from h = 1.5 μm to h = 14 μm. Though h = 14 μm is larger than a typical cell height, ridges of this height do not completely suppress crossing. We see that in these high ridges, the cell monolayer tends to become slightly undulated and threedimensional, with not all cells in contact with the substrate (SI Appendix, Fig. S10). Effects of this three-dimensional structure will not be fully captured in our 2D simulations. Increased ridge height of h = 14 μm leads to increased relative cell density near the +1 defect (Fig. 6 D and F ), consistent with our simulation results. Near the -1 defect densities remain low, with no clear systematic dependence on ridge height-similar to our simulation results for k r = 80-120. We see that by constraining cell migration and reducing ridge crossing, we increase density effects at the +1 defect core, as seen in our simulations, and as expected if cell division is driving the increase of density.

We find that 3T6 fibroblast alignment, velocity patterns, and density variation in our experiment are consistent with a model using shape-dependent cell division. We experimentally induced defects using ridges, which resulted in high fibroblast density near +1 defects and low density at -1 defects. However, unlike experiments with other cell types where cell accumulation at positive defects and depletion at negative defects were driven by collective migration (9-11), 3T6 fibroblasts did not manifest collective inward flow with highly aligned velocities. Instead, they moved in random azimuthal directions relative to the center of the +1 defect. This is consistent with earlier work arguing that fibroblast monolayers are less driven by migration and activity (13). To understand why density differences arise, we modeled fibroblasts as deformable elliptical cells. Our simulations found patterns of cell migration and alignment with the ridges that are consistent with experimental observations. Based on prior experiments on dependence of cell cycle progression on shape, we proposed a proliferation procedure where larger and more isotropic cells have higher probability to divide. This mechanism leads to density variations consistent with those in experiments. We predicted, using our simulations, that restricting cell movement across ridges would increase density at the +1 core-and confirmed it in experiments by modifying ridge height. Despite strong migration constraints, the marked density differences at defects were still present, which implies that cell division is important for explaining accumulation of cells.

Our model argues that the key factor controlling whether we get high or low density at a particular point is not the topological charge itself but whether cells can pack more efficiently isotropically or in an elongated state. In this sense, the specific geometry of confinement near the defect is important in determining the density. Supporting this idea, if we change the size of the -1 defect, we can develop a region of near-average density near the center of the -1 defect (SI Appendix, Fig. S11), though we still see overall that the -1 defect has a lower-thanaverage density. This hints that growth could be controlled both by shape and size of patterns. Further research is necessary to understand how different pattern variations impact growth.

Our model assumption that cell shape and size regulate division is consistent with past experiments (24,25). However, other experiments have argued that stress or pressure control proliferation and growth (26)(27)(28)(29)(30). These may be elements of a single-core mechanism, as cell shape, size, and stress are all intertwined (31).

Our results imply that patterned substrates can regulate the development of tissues via control of proliferation and not merely through controlling migration (32)(33)(34)(35). Similar approaches may help use mechanical cues to organize cells with limited motility. Beyond simply growth, other work shows that confinement and topology can provide cues to drive differentiation of cells (36,37). As cell area and aspect ratio are known to be important in determining the fate of individual cells (38)(39)(40), our work suggests that ridge patterns could be harnessed for controlled development (41, 42)-but the observed feedback between growth and cell shape means that computational modeling will be required to understand the effect of any given pattern. Changes of cell shape and aspect ratio are also seen in many patterning processes in development, including avian skin morphogenesis (43)(44)(45); control of division by local cell shape may allow for additional feedback between tissue growth and local alignment. Our results suggest that capturing the interplay of division, liquid-crystal alignment, and cell shape together is required to understand many patterning processes in eukaryotic development.

A. Simulation. We model cells as self-propelled elliptical particles with area A = π ab and aspect ratio AR = a/b where a and b are major and minor axis radii, respectively. The a and b can vary from cell to cell and will change over time; in our convention, a is chosen such that it is the larger axis of the cell, a ≥ b. We perform Monte Carlo simulations using the Metropolis method. Briefly, we propose changes to cell properties-these changes are accepted with probability min(1, e -E/T ) where E is the change in energy due to the proposed move and T is temperature. As in, e.g., the Cellular Potts Model and related models (46), this is not a physical temperature, but a value setting the likelihood of fluctuations of different sorts. We choose the temperature T = 1 to set the energy scale of the problem. The three central elements of the simulation are proposed moves, associated energies of the move, and cell division. A.1. Proposed moves. In one Monte-Carlo step (MCS), we iterate over all cells in random order and propose a single move for each of them. For a cell i, one of the four possible moves is attempted:

Rotate by φ i = 2δφ(X -0.5). 3. Change major axis radius by a i = δa(X -0.5).

where X, Y ∼ U [0, 1] are random variables sampled from the uniform distribution defined in unit interval [0, 1]. The parameters δr, δφ, δa, and δb represent the maximum possible displacement, rotation angle, change in semimajor axis length and change in semiminor axis length at each attempt, respectively (numeric values are given in Table 1). We reject or accept a move after each attempt based on energy change E i of the cell i due to the proposed move.

To calibrate timescales of cell growth and speed, each of the proposed move types has a different probability to be selected. We propose displacement with probability of 10%, rotation with 20%, and the two axis length changes each have 35% probability to be selected as a move. We have chosen this in part to ensure that cells quickly reach their steady-state shape, reflecting observations in experiments that, e.g., equilibration of fibroblast shape after division is much faster than significant motility (47). A.2. Cell energies. We accept moves following the Metropolis criterion, which depends on the change in energy from a move. The total energy of our system is composed of four distinct parts: geometric energies, cell-cell interaction energy, cell-ridge interaction energy, and motility energy. We describe each of these here.

Geometric energies penalize deviations from preferred size and shape. Cells have preferred area A pref and preferred aspect ratio AR pref . We penalize deviations from preferred values with energy cost quartic in relative deviations δ A = (A -A pref )/A pref and δ AR = (AR -AR pref )/AR pref :

Table 1.

1.2 μm Max. displacement proposal 2 deg Max. rotation proposal a 0.6 μm Max. semi-major axis change proposal b 0.6 μm Max. semi-minor axis change proposal A pref 1,400 μm 2 Preferred area AR pref 4 Preferred aspect ratio k A 18 T Area penalty strength k AR 4 T AR penalty strength k core 1,600 Tμm 4 Core penalty strength 0 0.5 T Gay-Berne interaction strength r maj c 100 μm Semi-major axis length of interaction cutoff r min c 43.75 μm Semi-minor axis length of interaction cutoff k move 1 T/μm Motility strength pol 30 min Polarity decay timescale flip 150 min Mean time between flips 1 Energy exponent 1 Energy exponent

where k A and k AR are area and aspect ratio penalty strengths. The shapes of energy curves are shown in SI Appendix, Fig. S12 A and B. Our goal in choosing these functions is to allow cells to easily change area and aspect ratio over a range of values close to their preferred values without significant energy cost. This reflects, e.g., for the area, that the cell can expand its height above the substrate, allowing it to make small area changes relatively easily. However, larger deviations result in substantial energy cost. The finite size of organelles and high nucleus stiffness relative to the cytosol implies that cells cannot be squeezed indefinitely. We model this feature via a core energy that introduces high energy cost if cell gets tiny but is much smaller when cells have a typical size,

whereaandbarethemajorandminoraxisradii.WeplotthiscurveinSIAppendix, Fig. S12C. Cell-cell interaction energy favors parallel alignment of long axes of cells. Cells interact with neighbors within cutoff region via a modified Gay-Berne potential that is extensively used in liquid crystal simulations. The potential favors mutual alignment of cells and it is strongly repulsive when cells are too close and weakly attractive if cells are separated by longer distances.

Here, we provide brief overview of the potential that we adapted for our simulations, detailed information can be found in refs. 15 and 48. The interaction depends on relative orientations of cells. We characterize orientation of a cell i at position r i = (x, y) by unit vector ûi = (u x,i , u y,i ) = (cos φ i , sin φ i ), where φ i is the angle the major axis of the cell makes with the x-axis of the simulation box. The interaction energy of a pair of cells located at positions r 1 and r 2 is given by

where û1 , û2 are cell orientations and r 12 = r 1r 2 is a vectorial distance between centers of cells (SI Appendix, Fig. S13A). The function r(û 1 , û2 , r 12 ) is a scaled and shifted distance between cells:

where σ 0 = 2b 2 1 + 2b 2 2 and σ (û 1 , û2 , r12 ) = 1/ r12 γ -1 r12 is an anisotropic range parameter that depends on size and orientations of cells via a matrix γ that depends on the size and orientation of both cells 1 and 2 as:

where l i = √ 2a i and d i = √ 2b i and the I is an identity matrix, and ûû indicates the dyadic product.

The term ε(û 1 , û2 , r12 ) = ε 0 ε ν a ε µ b is an anisotropic interaction strength. Here, ε 0 sets the general strength of interaction while ε a = 1/ 1χ 2 û1 • û2 and ε b = 1/σ (û 1 , û2 , r12 ) 2 scale strength based on size and relative orientation of cells. ν and µ are adjustable exponents set to 1 in our simulations and χ is given by

While the Gay-Berne potential of Eq. 5 has a long-distance attraction, we cut it off after a characteristic distance, reflecting that we do not expect cells to interact too far beyond contact. We compute the interaction energy of cell i with cells whose center is located within an elliptic area surrounding i. The elliptic area has same orientation as cell i and has semimajor and semiminor axis lengths r maj c and r min c , respectively (SI Appendix, Fig. S13A). This cutoff also allows us to speed up our simulations, as we only need to compute pairwise interaction energy between cells that are within this distance, which we track with a neighbor list. We update neighbor lists for all cells every time a cell divides or if any one of the cells moves by more than 25 μm with respect to its location during previous neighbor list update.

Ridge-cell interaction energy. Ridges are elevated with respect to rest of the substrate, so we penalize cell-ridge overlap. The energy cost of overlap is equal to the product of the ridge strength k r and fraction of the cell intersecting with the ridge, which we call ϒ:

[10]

To estimate the fraction of cell overlap ϒ we compute the overlap between individual points within the cell, where these points are sitting on three ellipses with same orientation and shape of the cell. The axis lengths of the outermost ellipse match the cell size (a, b), and the inner two ellipses have axis radii (2a/3, 2b/3) and (a/3, b/3), respectively. If the number of points on these three ellipses that overlap with ridges is N o and total number of points is N t , then the fraction is given by ϒ = N o /N t . Each "feeler" ellipse has 64 points separated evenly in polar angle (SI Appendix, Fig. S13C).

Cell motility energy promotes movements along long axis in the direction of polarity. Cells are animate entities that constantly convert chemical energy into mechanical movement. They often have persistent direction of motion that may change by itself or due to cues like electric field, chemical gradient etc. (19). The direction in which cell wants to travel is called (migrational) polarity, it points from rear of a cell where myosin contractions pull the "back" of the cell to the "front" where filopodia or lamellipodia push the frontier of the cell (49). We denote the polarity vector of a cell by p = (p x , p y ). In our model, when a cell rotates by φ, we correspondingly rotate the polarity. Because fibroblasts tend to move along the long axis of the cell (4), we choose the energy to promote motion along the long axis in the direction of polarity. For instance, if r = ( x, y) is a proposed displacement of one cell with polarity p then the motility energy change that results from this move is

where = (û • p)û is projection of polarity onto the long axis of the cell. Note that, in our approach, the magnitude of the polarity is irrelevant-only its direction contributes to the energy. k move sets the relevance of the motility compared to other driving forces. The energy we use here is akin to, e.g., energies used in the Cellular Potts Model to represent cell polarity and motility (23,50).

The polarity has positive feedback with the displacement in one MCS, r, which could also be zero if no displacement is proposed or the proposed displacement is rejected. At every MCS t, we update the polarity for each cell as in Szabo et al. (23):

where p t and p t-1 denote the new polarity at step t and polarity at previous step t -1, respectively, and τ pol is a polarity decay timescale parameter measured in units of MCS. If there were no displacement r = 0, polarity would decay in magnitude. However, because displacements tend to correlate with polarity and correlate more with polarity if many cells are locally pushing in the same direction, Eq. 12 tends to cause cells to locally align. The effect of changing polarity decay time has been systematically investigated in ref. 21. In this sort of model, increasing the time required for polarity to decay ensures that the polarity is largely controlled by the sum of previous displacements over a long time-generally making the migration more coherent.

In addition, since fibroblasts periodically reverse direction of motion (4), we stochastically flip cell polarity, p → -p. Reversal happens with probability 0.01 every 1.5 min (100 MCS). The number of tries needed for flip event has geometric distribution with success probability p = 0.01. Expected number of attempts needed for reversal is then 1/p = 100 (i.e., average flip time is τ flip = 1.5 × 100 = 150 min). This random flipping disrupts polar coherence, preventing cells from forming a uniformly rotating state as can be seen in, e.g., experiments on epithelial trains (14). A.3. Cell division. As in the monolayer experiments, in simulations our cells divide and proliferate. We initialize our system by putting circular cells of initial radius r 0 = 10 μm at random positions (excluding configurations with cellcell overlap) at a density of ρ ∼ 70 cells/mm 2 . We let cells evolve for 10 h without division, to ensure that they can relax to reasonable shapes. Then, we divide one cell every 1.5 min (100 MCS), choosing this rate to roughly match the experimental growth curve (SI Appendix, Fig. S2), and halt division once cells reach their terminal density of ρ f ∼ 2,000 cells/mm 2 (3,000 cells in our simulation box size of 1,200 μm × 1,200 μm). In our model, the number of cells as a function of time is always the same from simulation to simulation, but which cell divides at any point is stochastic. Cells have shape-dependent probability to be selected to divide: Given N cells in the simulation, cell i is selected to divide with probability

where α is a parameter that tunes how sensitive division probability is to cell shape and Z is a normalization factor chosen such that N i=0 p i = 1. When cell i of size (a i , b i ) and orientation ûi divides the two daughter cells, both will have size (a 1 , b 1 ) = (a 2 , b 2 ) = (0.4a i , 0.4b i ) and orientation û1 = û2 = ûi (SI Appendix, Fig. S13B). Divisions occur along the long axis of the cell; the choice of daughter cell sizes ensures that the aspect ratio of the cell is preserved. The choice of daughter cell size does mean that area is not conserved in the division-this is in part to avoid potential numerical problems with extreme cell-cell overlap which can be caused by division. We find that cells quickly grow up to a size comparable to nearby cells post-division when possible. A.4. Randomness and seeds. When varying parameters (e.g., in Figs. 4 and 6), we perform 100 simulations (indexed 1, 2, ..., 100) for each parameter set. To better understand the effect of the changed parameter, we keep the random number generator seed of each simulation fixed-so there are 100 distinct seeds, and when a parameter is varied, e.g., comparing α = 2 to α = 0 in Fig. 4, we are comparing simulations that have the same initial conditions and seeds. This choice is made to make it clearer that changes are systematically due to the effect of the parameter change and not randomness. However, because the set of 100 initial conditions are the same for all of our runs, we need to be confident that these initial conditions are not significantly driving our results. We provide a test of this in SI Appendix, Fig. S14, swapping initial conditions between the +1 and -1 defects. While there is some quantitative difference from earlier results in Fig. 4, we still see the key results that density is increased at the +1 defect as we make α → 2.

A.5. Broader modeling considerations. One contribution of our work is a cellbased framework, which can describe collective migration of highly anisotropic cells while resolving individual cell shapes and positions. This is one of many possible approaches to modeling collective cell migration (51,52). We discuss some of the broader choices we made here. We argue continuum tissue/active liquid crystal models (52,53) would be inappropriate to model these experiments, as continuum models are restricted to length scales that can average over many individual elements-incompatible with studying cells of typical width ∼20 μm in ridges with spacing 60 μm. Our approach also avoids issues with orientational anisotropy associated with lattice models like the Cellular Potts Model (46). In principle, phase field approaches (22,(54)(55)(56)(57)(58) would also avoid lattice artifacts, but our scale of ∼3,000 cells is an order of magnitude larger than typical applications of even simplified phase-field models (59)(60)(61). Earlier papers have modeled elongated self-propelled objects with particle-field and/or Gay-Berne approaches (62,63), though without explicitly describing deformability. Our approach is probably closest to the deformable self-propelled particle approach of Menzel and Ohta (64) and the deformable ellipsoids of Palsson and Othmer (65).

We have neglected in our paper the possibility that cells may create alignment of fibronectin or other extracellular matrix proteins that may play a role in longrange guidance (66). We were able to recapitulate alignment to ridges without this effect. However, it may be essential to understand longer-scale perfect alignment on unpatterned substrates (4). We also neglect potential couplings between cell shape and polarity (67-69), which can drive complex behaviors like cell circling and oscillation in response to fields, as observed in keratocytes (68). We have neglected these factors because we have no evidence that fibroblasts show circular migration behavior similar to keratocytes.

We have focused on how our division rules alter local effects of density in response to patterning. These division rules will also likely affect the mechanical properties of the tissue and the degree of fluidity (70,71). These are factors that might be important to study further in extensions to more motile cells than our 3T6 fibroblasts. Pressure feedbacks on growth rate and their interplay with fluidity have also been previously studied (72)(73)(74)(75).

B.1. Cell culture. The 3T6 mouse fibroblasts (from ATCC) are cultured in CellTreat tissue culture dishes using 90% Dulbecco's Modified Eagle's Medium (DMEM) [+] 4.5g/L glucose, L-glutamine, sodium pyruvate (Corning CellGro), and 10% Fetal Bovine Serum (Corning CellGro). When outside the incubator for longduration (>30 min) imaging, the growth medium is replaced with 90% CO 2 Independent Medium (Gibco) and 10% Fetal Bovine Serum (Corning CellGro), with 4.5 g/L L-glutamine added (Quality Biological). Cells are utilized only up to generation 20.

Growth curve. The experimental growth curve of SI Appendix, Fig. S2A is obtained using a standard method of seeding cells onto five Petri dishes with a density of 60 cells/mm 2 . To measure this seeding density, a subset of the suspended cells are stained using Trypan Blue and counted using a hemacytometer (10 μL volume). Every day, one dish is selected, for which cells are resuspended and counted with the same method. Cells in the five dishes are counted after 24, 48, 75, 95, and 120 h, respectively. B.2. Substrate preparation. The topographic features are patterned using a mold of SU-8 on glass. SU-8 is a negative photoresist, a hard polymer which crosslinks by exposure to UV light. To create 1.5 μm-tall ridges, we use SU-8 2002 (MicroChem), by spincoating the SU-8 at a maximum speed of 3,000 rpm for 30 s. To create 6 μm-tall ridges, we use SU-8 2005 (MicroChem), by spincoating the SU-8 at a maximum speed of 3,000 rpm for 30 s. To create 14 μm-tall ridges, we use SU-8 2010 (MicroChem). In this case, we use a maximum speed of 2,000 rpm for 60 s.

Polydimethylsiloxane(PDMS)substratesarepreparedfromSylgard184(Dow Corning) with 15% curing agent. After mixing, the PDMS is degassed at room temperature and then poured over the SU-8 patterned substrates. These are cured in a vacuum oven at 60 • C for 2 h. The patterned PDMS is then prepared for cell culture.

We submerge the patterned PDMS in ethanol for 20 min to sterilize it. The substrates are then dried at room temperature, and then coated with a minimal volume of fibronectin from bovine plasma (MilliporeSigma, 25 μg/mL in PBS) for 45 min at room temperature prior to use for cell culture.

Prior to seeding onto the patterns, a subset of the suspended cells are stained using Trypan Blue and counted using a hemacytometer (10µL volume). Cells are seeded at a concentration of 5 × 10 5 cells/dish, with each dish being circular with a diameter of 100 mm. This corresponds to 60 cells/mm 2 .

We note here that the data for cells growing on 1.5 μm shown in Figs. 1D and 6F are from our previous publication (12), For these datasets, the seeding density was not held constant.

As a control, the ridge height and pattern quality of the PDMS with ridges up to 6 μm are verified using a Keyence VK-X200K color 3D laser scanning microscope. The 14-μm ridge heights are verified by imaging the cross-section of the PDMS using a 50X objective on a Nikon LV Pol microscope with a Nikon DS-Ri2 camera. For all ridge heights, we analyze the images in the acquisition software by measuring the height of the top surface at a series of locations both on and off the ridges and computing the average difference, with SE as uncertainty. Measurements of ridge height are shown in SI Appendix, Fig. S15.

Preparation/Staining. The cells' nuclei are stained using Hoechst 33342 dye (10 mg/mL stock solution, Invitrogen), diluted at a ratio of 1:1000, followed by 15 min in the incubator. When they are stained for the purpose of acquiring a video, the dye solution is removed, and the dish is filled with CO 2 -independent media.

For the samples with 1.5 μm high ridges, there are three images each of cells around +1 defects and -1 defect of which the nuclei are stained using NucRed Live 647 (Invitrogen) rather than Hoechst 33342. The dye is added following the suggested protocol of 2 drops/mL of media, followed by 15 min of incubation. Data about nuclear orientation and cell density from these images are incorporated into Figs. 1 C and D, Fig. 3B, and Fig. 6F.

For confocal images, cells are fixed before being stained with Hoechst 33342 and Phalloidin-iFluor 594 conjugate (AAT Bioquest). For fixing, the cells are first incubated in 4% paraformaldehyde for 10 min at room temperature. Then, they are washed with PBS three times, each time incubated at room temperature for 5 min. Then, they are incubated in PBS with 0.1% Tween-20 at room temperature for 10 min, followed again by washing three times in PBS for 5 min each. Once the cells are fixed, they are stained with a solution with 10 μg/mL Hoechst 33342 and a working solution of the Phalloidin stain. This Phalloidin stain is composed of 1 μL of Phalloidin-iFluor 594 Conjugate solution (AAT Bioquest) diluted in 1 mL of PBS with 1% Bovine Serum Albumin. The cells are incubated in this solution for 20 min prior to imaging.

Microscopy. Phase contrast and fluorescent imaging in 2D is done using a Nikon TI-Eclipse microscope using a Hamamatsu Orca-flash camera. Large format images are acquired by translating the stage with 15% overlap between adjacent frames, and then, the images are stitched together using Stitching (Grid/Collection stitching) plugin in ImageJ (76). At each location, we take one phase contrast and one fluorescent image.

For video acquisition, phase contrast and images of nuclear fluorescence are taken at 6 min intervals. While acquiring the video, the Tokai Hit ThermoPlate microscope stage is heated to 37 • C, and the stage is covered with a plastic sheet. Twice a day, or as needed, fresh CO 2 -independent media is added to refill the dish, which loses medium due to evaporation. Autofocus is performed in NIS-Elements (Version 5.02.01) before each image is captured, and the light is switched off between frames.

Confocal microscopy. Toquantifytheprevalenceofcellsgrowingoverridges, the cells are imaged using a Leica SP8 confocal microscope with a White Light Laser and Leica HyD detectors and 40X objective. They are imaged after being fixed, permeabilized, and stained with Hoechst 33342 dye (10 mg/mL stock solution purchased from Invitrogen), diluted at a ratio of 1:1,000. For each step, we take a stack of images at different heights from three channels: a bright field channel, a channel collecting the information from the nuclear fluorescence, and a channel collecting the information from the actin filaments stained with Phalloidin 594.

Confocal microscopy. To measure the fraction of nuclei that are on or off the ridge f on , f off as used in Fig. 6 D and E, images are segmented with the software IMARIS 9.8.2 using the "Create Surface" feature to identify the cells growing over the ridge from those growing between ridges. In the bright field images, we identify clearly defined lines on the bottom plane of the z-stack corresponding to the location of the edges of ridges. Using a drawing tool in IMARIS, we generate a shape on each edge. Then on the top plane of the z-stack, we paste the same shapes that were generated on the bottom surface. Around the circular ridges used to generate +1 defects, this demarcates a cylindrical shell which we identify as the "on wall" region of the image. This procedure is repeated for every ridge in every image. Once these regions are demarcated, we use the command "Mask Selection" to create a new color channel containing only voxels from the fluorescent image which are "on wall." These regions, identified with the method described above, are used to create a mask on the corresponding images of fluorescent nuclei.

Once the masks are made, the 3D images are analyzed using ImageJ. First, all the channels are analyzed separately; then, the images are collapsed to 2D using the maximal intensity of a voxel at a given xy-position. An intensity threshold is applied to create binary images, only showing the nuclei at given xy-positions. Fractions of black/white pixels on/off walls are computed to identify the area fractions on/off walls occupied by nuclei. From these, we report the average area fractions and the SDs.

Cell alignment. Quantification of cell orientation also followed the protocol of ref. 12. Briefly, we define the axis along which the topographic features orient the cells as a function of its azimuthal coordinate with respect to the defect, following the protocol described in ref. 12. For a +1 defect, this is 90 • more than the coordinate itself, modulo 180 • . For a -1 defect, this is 180 • minus the coordinate in the first and second quadrants, and 360 minus the coordinate in the third and fourth quadrants. We then compare the deviation in the direction of the cell's major axis from this axis, which is the deviation in the cell's alignment from the expected or patterned angle. The orientation of the cell's major axis is determined by fitting ellipses to nuclei identified in the fluorescent images in ImageJ, first by creating a binary mask of the image.

Cell density vs. distance from defect. We identify the center of the nearest patterned defect from the phase contrast images. For each pixel in the fluorescence image of the nuclei, we compute the distance from the nearest defect. We then compute the sum of the intensities and the number of pixels in each 60 μm ring out to a distance of 600 μm, and from these values, we compute an average nuclear intensity in each ring. In each plot, we report the average and the SE, after normalizing by dividing by the average intensity measured within 600 μm from the defect center.

Cell tracking. Cell tracking is achieved using the TrackMate plugin in ImageJ (77). The simple Linear Assignment Problem (LAP) tracker is used, which does not detect merging and splitting events. The trajectories of each cell are imported into Matlab for further analysis. Fig. 1G shows the direction of displacement of the cells in the first hour of the video. Only the cell paths that are continuously identified for every frame in the first hour are included in the image. This includes the majority but not all of the cells in the frame. Then, these cells are identified as moving clockwise or counterclockwise, and a vector with a fixed length is drawn on the first frame of the video. Its direction indicates the direction (but not the magnitude) of the net displacement of the cell in the first hour and its color indicates whether the cell moves clockwise or counterclockwise.