A predictive model of intercellular tension and cell-matrix mechanical interactions in a multicellular geometry

Epithelial cells form continuous sheets of cells that exist in tensional homeostasis. Homeostasis is maintained through cell-to-cell adhesions that distribute tension and balance forces between cells and their underlying matrix. Disruption of tensional homeostasis can lead to Epithelial-Mesenchymal Transition (EMT), which is a transdifferentiation process in which epithelial cells adopt a mesenchymal phenotype, where cell-cell adhesion is lost and individual cell migration is acquired. This process is critical during embryogenesis and wound healing, but is also dysregulated in many disease states. To further understand the role of intercellular tension in spatial patterning of epithelial cell monolayers, we developed a multicellular computational model of cell-cell and cell-substrate forces. This work builds on a hybrid Cellular Potts-finite element model to evaluate cell-matrix mechanical feedback of an adherent multicellular cluster. Thermodynamically-constrained cells migrate by generating traction forces on a finite element substrate to minimize the total energy of the system. Junctional forces at cell-cell contacts balance these traction forces, thereby producing a mechanically stable epithelial monolayer. Simulations were compared to in vitro experiments using fluorescence-based junction force sensors in clusters of cells undergoing EMT. Results indicate that the multicellular CPM model can reproduce many aspects of EMT, including epithelial monolayer formation dynamics, changes in cell geometry, and spatial patterning of cell geometry and cell-cell forces in an epithelial colony. Author summary Epithelial cells line all organs of the human body and act as a protective barrier by forming a continuous sheet. These cells exert force on both their neighboring cells as well as the underlying extracellular matrix, which is a network of proteins that creates the structure of tissues. Here we develop a model that encompasses both cell-cell forces and cell-matrix forces in an epithelial cell sheet. The model accounts for cell migration and proliferation, and regulates how cell-cell adhesions are formed. We demonstrate how the interplay between cell-cell forces and cell-matrix forces can regulate the formation of the epithelial cell sheet, the organization of cells within the sheet, and the pattern of cell geometries and cell forces within the sheet. We compare computational results with experiments in which epithelial cell sheets are disrupted and cell-cell junction forces are measured, and demonstrate that the model captures many aspects of epithelial cell dynamics observed experimentally.

The epithelium is characterized by polarized sheets of cells that form by 2 self-organization and reside in a mechanical equilibrium (reviewed in [1]). This 3 mechanical equilibrium is maintained by regulation of both adhesion between 4 neighboring epithelial cells (cell-cell) as well as adhesion between epithelial cells and the 5 underlying extracellular matrix (cell-matrix). Cells generate cytoskeletal tension via 6 actomyosin contractility, which is transmitted to the underlying matrix, while cell-cell 7 adhesion mechanically couples abutted cells and distributes cytoskeletal tension to 8 neighboring cells. This physical cellular interconnectivity and balance of tension at the 9 cell-matrix and cell-cell interfaces produces a coupled monolayer that acts as a cohesive 10 structure in static equilibrium. 11 Maintenance of static equilibrium in the epithelial sheet is essential to maintaining 12 barrier and signaling functions of the epithelial sheet; however, disruption of the static 13 equilibrium plays an important role in both physiological phenomena such as 14 embryogenesis and pathological states including fibrosis and tumorigenesis [2,3]. 15 Mechanical equilibrium relies on tissue scale coordination of mechanical dynamics 16 extending beyond local cell-cell and cell-matrix adhesions [4]. Local perturbations to the 17 equilibrium state result in localized tension in the monolayer and a disruption to the 18 equilibrium. For example, the cellular phenomena known as epithelial-mesenchymal 19 transition (EMT), which is essential for embryogenesis and tissue morphogenesis, but 20 which has also been implicated in tumorigenesis and fibrotic diseases, is initialized by 21 perturbations in cell-cell adhesion. This process results in a phenotypic switch in which 22 epithelial cells transdifferentiate into mesenchymal cells (reviewed in [5]). The 23 perturbation in cell-cell adhesion redistributes tension in the monolayer, and cell-matrix 24 adhesion compensates for the resulting localized stress. As such, spatial patterning of 25 mechanical stress can facilitate phenotypic regulation and is crucial to both 26 maintenance and disruption of tissue homeostasis [2,4,6]. 27 Previous studies have explored the role of cell-cell adhesion in maintaining tensional 28 homeostasis in the epithelial monolayer: increasing cellular contractility has been shown 29 to stimulate formation of cell-cell junctions [7], and subsequent transfer of force to the 30 cell-cell adhesion allows for stress distribution about the monolayer to maintain 31 tensional homeostasis [4,6]. As a result, mechanical gradients form that define spatial 32 patterns and provide positional information within the monolayer. Both in vitro and in 33 silico studies have demonstrated that the forces of a monolayer correspond to its 34 geometry [8,9]. 35 In this work, we explore the role of cellular adhesion in maintaining tensional 36 homeostasis of epithelial monolayers. To simulate epithelial monolayers, we extended a 37 model developed by van Oers et al, which consists of a hybrid Cellular Potts model 38 (CPM) and finite element model (FEM) [10]. The model simulates individual cellular 39 traction forces based on their geometric size and shape, as has previously been modeled 40 and validated by one of the senior authors of this work [11]: cellular traction forces are 41 proportional to the first moment of area (FMA) about each point in the individual impose a thermodynamic constraint and govern the dynamics of individual cells in the 46 CPM. In the current work, we incorporate the formation of cell-cell adhesion between 47 neighboring cells to accurately represent the biology of epithelial cells. We extend the 48 Lemmon and Romer FMA model to multicellular clusters, and model traction forces 49 based on the multicellular geometry rather than the individual cell. Thus, individual 50 cell traction forces are proportional in magnitude to the distance from the centroid of 51 the multicellular cluster, instead of the centroid of the individual cell. 52 In the original Lemmon and Romer model, each cell is in static equilibrium: because 53 traction forces are proportional to the first moment of area, and the centroid by 54 definition is the point where the integral of the first moment of area is zero, all traction 55 forces within a cell must sum to zero. However, when we calculate traction forces based 56 on the multicellular cluster, each individual cell is no longer in static equilibrium. 57 Previous studies have suggested that cells in epithelial monolayers exist in a 58 quasi-equilibrium, even when cell-cell junctional forces are present [7]. As such, we 59 model the force applied to the cell-cell junction as the balancing force that opposes the 60 traction forces for that cell, resulting in a quasi-equilibrium for each cell. This 61 assumption has been observed experimentally in epithelial cell pairs, in which the 62 junction force is equal and opposite to the net traction force [7]. We thus are able to 63 predict the formation of an epithelial monolayer, including epithelial cell geometry,  Simulations demonstrate that traction forces of multicellular colonies scale linearly 76 with the size of the colony, independent of the individual cell geometry. Additionally, we 77 demonstrate that the model can be generalized to predict the distribution of junctional 78 forces across a monolayer: junction forces are predicted by a quadratic function that is 79 highest at the colony center and decays towards the cell edge. These predictions are 80 independent of indvidual cell geometry and are consistent with existing literature [12].  To expand this model to adherent cell monolayers, we incorporated several 88 advancements: first, cellular traction forces were predicted from the FMA model [11] 89 based on a cell cluster geometry, not on individual cells. As such, cells in contact with 90 neighboring cells "adhere" and begin to generate traction forces as a cohesive unit.

91
Second, we assume that each cell in a multicellular cluster still maintains a static 92 equilibrium, as has been suggested previously [7]. As such, we require the force acting 93 on cell-cell junctions to counter the net traction force for each cell, as illustrated in a 94 simple two cell example (Fig 1C, left).  strains (black vectors) for two scenarios. In the first, traction force is calculated from 97 the first moment of area (FMA) about the single cell geometry and each cell is in static 98 equilibrium. As a result, the net imbalance for each cell is zero and no force is 99 transferred across the cell-cell junction (Fig 1A). In the second scenario, traction force is 100 calculated from FMA about the multicellular geometry and each cluster is in static 101 equilibrium ( Fig 1B). The net force imbalance for each cell is balanced by the 102 intercellular tension, which transfers the traction force to neighboring cells. Without 103 redistribution of cytoskeletal stress to neighboring cells across cell-cell junctions, cellular 104 alignment is localized and multicellular structures behave as partially cooperative 105 networks with discordant substrate strains ( Fig 1A, S1 Video), as demonstrated by van 106 Oers et al [10]. In contrast, traction force distribution across cell-cell junctions to 107 neighboring cells results in highly cooperative networks with a uniform spatial gradient 108 of substrate strains. The formation of these cohesive multicellular clusters resembles an 109 epithelial monolayer with preferential localization towards the boundary (Fig 1B, S2 110 Video). In the resulting multicellular clusters, net traction forces have a magnitude and 111 direction at any given point proportional to the FMA about that point in the cluster,  With the key addition that traction forces are governed by the FMA model about the 135 cluster geometry rather than the single cell geometry, the previous results illustrate 136 distinct spatial patterning representative of epithelial monolayers. We next utilized our 137 model to simulate epithelial monolayer and associated EMT-like dynamics. One key 138 aspect of the epithelial phenotype is contact inhibition: that is, the propensity of a cell 139 to stop migration when a neighboring cell is encountered [14,15]. As epithelial cells 140 undergo EMT and become more mesenchymal, contact inhibition is reduced [16]. To i.e., a higher J cc /J cm ratio reflects increased contact inhibition between adjacent cells. 148 For each simulation, we measured the steady-state monolayer confluence, average cell 149 area, total cell count, and relative net cellular traction forces, averaged over 5 150 simulations with distinct random cell seeding, and plotted these measures as a function 151 of the J cc /J cm ratio. These simulations were then repeated for 3 distinct values of Results indicate similar trends between the single cell and multicellular FMA models, 154 with the exception of net cellular traction force, which must be zero for a cell in static 155 equilibrium in the single cell FMA model (Fig 3D). Beyond a critical point 156 (J cc /J cm = 2), high cell contact inhibition precludes the formation of confluent 157 monolayers ( Fig 3A, E). Further, we find that the time course of monolayer confluence 158 only weakly depends on cell contact inhibition below this critical point, i.e. for 159 conditions that form confluent monolayers ( Fig 2C). Similarly, increasing cell contact 160 inhibition results in smaller cell area (Fig 3B, F) and higher cell count (Fig 3C, G). In 161 the multicellular FMA model, net traction force per cell decreases as the J cc /J cm ratio 162 increases. We find that higher substrate inhibition, i.e., increased J cm , tends to increase 163 the sensitivity to the J cc /J cm ratio for all measures. Thus, these data indicate that a 164 loss of contact inhibition leads to larger cells, lower cell count, and in extreme cases, loss 165 of confluence. phenotype. Together, these results indicate that this parameter may serve as a suitable 175 comparison to in vitro models of growth factor induced EMT. We thus compared these 176 results to experiments in which EMT was induced by the soluble growth factor TGF-1, 177 as has previously been detailed [17]. Representative immunofluorescence images of 178 MCF10A cells treated with increasing dosages of TGF-1 illustrate a phenotypic switch 179 from cortical actin, which is typically observed in epithelial cells, to pronounced actin 180 stress fibers associated with the mesenchymal phenotype ( Fig 4A). In these confluent To examine spatial trends, we segmented the simulation domain into a 5 x 5 grid of bins 195 and calculated the mean junction force magnitude within each bin ( Fig 5E). The spatial 196 distribution of junction forces is pronounced, with the largest forces in the interior and 197 smallest in the corners ( Fig 5F). However, interestingly, we find minimal variation in the 198 spatial trends between low, medium, and high contact inhibition ratios. Simulated intercellular tension is depicted as the net magnitude for high, medium, and low interaction energy (J cc /J cm ) ratios. (E) Intercellular tension magnitudes are shown as a 5 x 5 grid with (F) their associated bar graphs averaged at the corners, edges, and interior; n=5, * with line denotes significance between each location.
We next sought to compare these with experimentally-measured junction forces. To 200 measure cell-cell junction forces experimentally, Madin-Darby Canine Kidney Cells 201 (MDCKII) cells were stably transfected with a full-length E-cadherin force sensor, as 202 previously described [18]. Briefly, the force sensor consists of two fluorophores coupled 203 by a polypeptide that exhibits elasticity. The two fluorphores are designed such that, 204 when in close proximity, the pair exhibits Forster Resonance Energy Transfer (FRET): 205 that is, emission light from the first fluorophore is absorbed by the second fluorophore, 206 which emits light. As the sensor is stretched and the fluorophore pair moves apart, the 207 excitiation of the second fluorophore by the first fluorophore decays, resulting in a loss 208 of FRET excitation relative to excitation of the first fluorophore. This force sensor was 209 inserted into E-Cadherin, which comprises the homophilic binding event in cell-cell 210 junctions known as adherens junctions. Validation and functionality of this sensor has 211 been previously demonstrated [19,20]. EMT was again induced by increasing dosage of 212 (TGF-1) (Fig 5A). FRET ratio reflects the energy transfer between the two 213 fluorophores, in which FRET ratio is inversely proportional to tension on the FRET uniform low FRET ratio, indicating high cell-cell tension throughout the monolayer (Fig 220  5B). TGF-1 treatment increased FRET ratio, indicating a drop in overall tension.

221
Additionally, a small spatial gradient was established, with higher FRET ratios (lower 222 cell-cell tension) in the corner and edges and lower FRET ratios (higher cell-cell tension) 223 in the interior of the monolayer, consistent with a spatial gradient of larger junction 224 forces in the center and decreasing towards the edges and corners (Fig 5C). 225 Thus, we find that simulated cell-cell junction forces predict a spatial trend of TGF-1 accentuates this trend, resulting in a large spatial gradient in cell area ( Fig 6B). 247 In contrast, simulated cell area exhibited substantially reduced spatial variation 248 compared to experimental cell area ( Fig 6C). Furthermore, the effects of contact 249 inhibition had a relatively minimal effect on spatial variation of cell area, resulting in 250 slightly increased cell area at the monolayer interior ( Fig 6D). Thus, the lack of 251 accounting for heterogeneous cellular properties, specifically cell area, is a key limitation 252 of our model. Since cells undergo profound phenotypic changes throughout EMT, it 253 would be reasonable that these changes lead to parameter changes within the CPM for 254 each individual cell; incorporating these changes in cell phenotype into the CPM 255 component is a primary future goal for the model development. such that ! J n 1,n = ( µLf (n 1 2 ), 0).

287
Next considering forces on cell n 1, the junction force from cell n 2 to cell n 1 288 must balance both the net traction force ! T n 1 = ( µLf ((n 1) 1 2 ), 0) and junction 289 force ! J n 1,n , such that ! J n 1,n 2 = (µLf (2n 2), 0). Similarly, junction force from 290 cell n 3 to cell n 2, ! J n 2,n 3 = (µLf (3n 9 2 ), 0). In general, we can show that the 291 intercellular tension from cell k to k + 1, Thus, the junction force at the center onto the left edge of cell 1, drop-off (due to the k 2 term in the magnitude of ! J k+1,k ) that is predicted as junction 296 position k increases towards the periphery. A representative example of the CPM model 297 illustrates the distribution of traction forces (blue) and junction forces (red) in a 298 confluent monolayer (Fig 7B) and both the linear increase in traction force magnitude 299 July 9, 2019 11/27 from the monolayer centroid and the quadratic drop-off in junction force magnitude 300 (Fig 7B).

301
Thus, for a monolayer of a given size, i.e., fixed T , Eqn. 1 predicts that for a smaller 302 cell size (decreased L and thus increased n), the magnitude of junction forces are larger 303 throughout the monolayer, which is consistent with experimental measurements of lower 304 FRET ratios (i.e., higher tension) in non-treated epithelial monolayers (Fig 5C). 305 Further, in TGF-1-treated monolayers, more mesenchymal-like larger cells at the 306 monolayer periphery would be expected to have more focal adhesions per cell, in 307 contrast with epithelial-like smaller cells in the interior. Additionally, while larger cells 308 at the periphery will reduce junction forces locally, due to the cumulative nature of 309 junction forces required to maintain mechanical equilibrium originating at the periphery, 310 this local reduction in junction forces would be expected to have a greater influence on 311 interior junction forces. All of these considerations would be predicted to reduce the 312 magnitude of the spatial gradient, also consistent with smaller spatial gradients observed 313 experimentally. Thus, we expect that our future work incorporating spatial variations in 314 cell size in the CPM model will more accurately reproduce experimental results. 315 We can further generalize this example and consider the continuous limit in the 316 spatial dimension, in which the traction forces ⌧ (x) in the x-direction at position x (for 317 x > 0) are given where (x) is the spatial distribution of focal adhesions per unit length. Junction forces 319 J(x) at position x are then by definition the second moment of area, evaluated from the 320 cluster periphery T to position x, where again x = 0 corresponds with the cluster center, 321 For uniform focal adhesion distribution, (x) = f /L, we can integrate Eqn. 2, and using 322 the relationship x = kL, the result is equivalent to Eqn. 1.

324
In this study, we illustrate a generalized framework for predicting the spatial tension redistributes from the cell-cell junctions to the cell-matrix attachments, which 347 allows for increased mobility, growth, and spreading [15]. Our model represents this shift 348 by altering contact penalties within the cell-cell and cell-matrix interaction energies. By 349 altering the cell-cell contact energy, the model captures the contact inhibition of 350 neighboring cells in vitro. However, simulating EMT via changes in the contact energy 351 is not sufficient to capture all dynamics: in the CPM model, a defined value for optimal 352 cell area constrains the simulated cell area that, in turn, limits cell-matrix adhesion.

353
The shift from cell-cell contact to cell-matrix adhesion is indirectly restricted as a result. 354 The spatial distribution of intercellular tension therefore predicts the spatial distribution 355 of cell area, which would seem to indicate a shift towards cell-matrix adhesion. Future 356 work will account for phenotype-dependent changes in the optimal cell area. dynamics such as morphogenesis [21][22][23]. Agent-based models have been utilized to 362 study cellular remodeling in response to mechanical perturbations, such as infarcts and 363 wound healing [24][25][26]. The CPM framework has also been utilized to study cell-matrix 364 interactions via extracellular matrix remodeling, in settings such as metastatic cancer 365 cell migration and angiogenesis [27][28][29].

366
Our work builds on prior studies from Merks and colleagues that have demonstrated 367 how local mechanical interactions can drive global cellular patterning and structure, 368 using a hybrid CPM-FEM framework [10,30,31]. Multiscale modeling studies from 369 Chaplain and colleagues have predicted that junction forces are redistributed as cells 370 form colonies, which in turn can drive intracellular signaling pathways [32][33][34]. information within a monolayer that regulates cellular phenomena, such as cell growth, 376 proliferation, and migration. This is of particular interest to spatial regulation of EMT, 377 during which cell stress is distributed to the monolayer periphery [35]. Connecting 378 biochemical and mechanical signaling, the dependence on E-cadherin further suggests 379 that intercellular tension may serve as a predictor of EMT.
and Boltzmann statistics determine the probability of a possible lattice configuration where H is the Hamiltonian defined in Eq. 4 and T > 0 is a temperature term that 412 captures intrinsic cell motility.

413
The area term H area approximates the cell area constraint as a deviation of the cell 414 area relative to the target area such that where a( ( ! x )) is the area of a given cell determined by number of lattice sites occupied 416 by that cell, A 0 = 312.50 µm 2 is the target area for all cells, and area = 500 is an 417 elasticity coefficient that maps deviations from the target area to a magnitude of energy. 418 The contact term H contact represents costs due to contact between neighboring 419 pixels, with different energies associated with cell-cell and cell-matrix interfaces: where J( ( ! x ), ( ! x 0 )) defines the interaction energy between adjacent lattice sites  Lastly, the durotaxis term H durotaxis introduced in van Oers [10] mimics the 425 tendency for cell migration along gradients of mechanical strain. In particular, this term 426 captures preferential cellular extension into lattice sites of higher strain July 9, 2019 14/27 The durotaxis = 1 term determines cell sensitivity to durotaxis; g( ! x , ! x 0 ) is 1 if a cell 428 extends into a target site ! x 0 and -1 if a cell retracts; and the v 1/2 · v m are defined such 429 that extension and retraction are greatest parallel to the major and minor principal 430 strain axes, v 1 and v 2 respectively, and negligible perpendicular to it. The sigmoid 431 function h (E) captures the preference for stiffer substrates, which assumes this preference has a minimal stiffness for spreading and reaches a 433 maximum ↵ = 10 at rate = 5 ⇥ 10 4 kPa 1 and the half-max stiffness as 434 E ✓ = 15 ⇥ 10 3 kPa. E(") is the cell perception of substrate strain stiffening, where " st = 0.1 determines the rate of strain-stiffening, " is the substrate strain, and To describe the substrate strain that governs durotaxis, we assume that a uniform, 441 isotropic, and linearly elastic two-dimensional substrate deforms to cellular traction 442 forces projected from the CPM (described below). The CPM lattice is mapped to the 443 finite element lattice by relating each CPM lattice element to a finite element node. We 444 solve the linear system for the displacement u at each node, where K is the global stiffness matrix assembled where B is the conventional strain-displacement matrix and D is the material property 450 matrix under plane stress conditions relating the Young's modulus, E = 10 kPa, and Poisson's ratio, v = 0.45, assuming 452 planar stress. Lastly, B relates the local node displacements to the local strains by 453 " = Bu n (14) in which " is a vector of the strain tensor ".

454
Traction forces 455 Prior work of van Oers and colleagues [10] assume that individual cell geometry relates 456 to traction forces in the CPM by the first moment of area (FMA). Application of the 457 FMA model to single cell geometries is previously described by one of the senior authors 458 of this work [11]. In brief, the single cell FMA model assumes that each node i in a 459 CPM cell exerts a force on all other nodes j in the same cell that is proportional to 460 the distance between those nodes where µ is a scaling factor that relates cell geometry to traction forces. For simplicity, 462 we assume µ = 1 nN µm 1 and report forces as relative arbitrary units (a.u.). As shown 463 in Lemmon and Romer [11], the resulting traction force at each CPM node is directed 464 towards the cell centroid with magnitude proportional to the distance from the node to 465 the centroid.

466
Here, we extend these previous works of the FMA model to describe the magnitude 467 and direction of traction forces acting about a point in a multicellular geometry. For the 468 multicellular FMA model, we assume that the boundary of two cells constitutes a 469 cell-cell adhesion such that two or more adjacent cells behave as a single structural unit 470 or cluster. We define an adjacency matrix A, where A is a N cell ⇥ N cell matrix, such cell T within a cluster may not be equal to 0. Adapting a recent approach by Ng and 486 colleagues [38], we hypothesize that junction forces are a reaction force, balancing the 487 net traction force to maintain static equilibrium of each cell in a multicellular cluster.

488
The multicellular FMA model is applied to calculate T for each cell, and then we 489 impose mechanical equilibrium on the multicellular clusters by relating the traction 490 force to force across the cell-cell adhesion, such that for all cells , where n defines the set of "neighbors" of cell , i.e., A , 0 = 1, and J , 0 is the 492 junction force from cell 0 to cell (see S2 Fig). Eq. 16 defines N cell linear equations, 493 with N 2 cell unknown J , 0 terms. We further constrain the junction force calculations by 494 assuming that junction force pairs are equal in magnitude and opposite in direction, i.e., 495 July 9, 2019 16/27 for all ( , 0 ) such that A( , 0 ) = 1.

496
Combining Eqs. 16 and 17, we arrive at a linear system with a set of N cell + N junc 497 equations and N 2 cell unknowns (see S2 Fig), where N junc is the number of intercellular 498 junctions, which can be determined by the sum of the terms above (or below) the main 499 diagonal of A, with a maximum value of N cell (N cell 1)/2. In practice, linear systems 500 for Eqs. 16 and 17 are determined separately to both the xand y-components of the 501 traction and junction forces.

502
For nearly all cluster arrangements, the resulting linear system is overdetermined.

503
Analogous to the CPM thermodynamic energy minimization, we assume that the 504 solution to be the minimization of junction force for each cell pair in the cluster, such 505 that J , 0 terms are calculated as the minimum norm least-squares solution to the linear 506 system (using the MATLAB lsqminnorm function).

507
Cell division 508 We incorporate cell division into the CPM model to reproduce epithelial cell capacity to 509 proliferate and form a confluent monolayer. For simplicity, we assume that if an substrate by identifying the corresponding nodes. At a given instant, the single cell or 523 multicellular geometry is sufficient to define cellular traction forces at each node, using 524 the single or multicellular FMA model, as described above, respectively. The resulting 525 traction forces govern the displacement at each node and determines the strain in the 526 finite element mesh, which in turn is used in evaluating H durotaxis .

527
Cell movement consists of copy attempts of randomly selected pixel at each Monte 528 Carlo step (MCS). For each pixel to have equal probability of selection, each MCS has a 529 total of 10 4 copy attempts. For each copy attempt, a pixel is selected and randomly 530 perturbed; the sum of interaction energies with each pixel in the Moore neighborhood, 531 P J( (x, x 0 )), determines the H contact term. Lastly, the cell area before and after the 532 copy attempt provides the H area term. Together, the net change in the Hamiltonian 533 associated with that copy attempt, i.e. H( (x, x 0 )), provides the local energy for the 534 cell before and after the copy attempt. The copy attempt is accepted ( (x) ! (x 0 )) 535 with probability determined by the partition function (Eq. 5) for H > 0 and 536 probability 1 for H < 0.

537
For parameter analysis, the parameter set consisted of each combination of cell-cell 538 interaction energies and cell-matrix interaction energies, J cc and J cm , respectively, each 539 repeated with a uniquely seeded random number. The confluence is determined by the 540 ratio of total cell occupied pixels to the total grid area. The cell area is number of pixels 541 occupied by each unique cell state, and the cell count is the number of unique states.   Microcontact printed square islands were generated as previously described [39]. Briefly, 562 250 µm x 250 µm squares were constructed by generating a negative mold template on 563 a silicon wafer made from an epoxy-type, near-UV photoresist (SU-8; Microchem) using 564 traditional photolithographic techniques. A replica-mold of poly(dimethylsiloxane)  Immunofluorescence microscopy 574 MCF10A and MDCKII cells were plated on microcontact-printed laminin islands at cell 575 densities that resulted in near-confluent monolayers. After 6 hours, samples were rinsed 576 in culture medium to remove non-adherent cells. Cells were cultured for 18 hr and were 577 then transferred to EGF-and serum-free culture conditions for 2 hr to induce an 578 epithelial phenotype. Cells were then incubated with or without TGF-1 for an 579 additional 48 hours. Cells were permeabilized with 0.5% Triton in 4% paraformaldehyde 580 for 2 minutes, then incubated in 4% paraformaldehyde for 20 minutes. Several 581 PBS-rinses were performed, followed by blocking in 0.1% BSA and labeling with primary 582 antibody for 30 minutes at 37 degree C. Cells were then blocked again in 0.1% BSA and 583 incubated with the appropriate secondary antibody for 30 minutes. Images were 584 acquired on a Zeiss AxioObserver Z1 fluorescence microscope using ZEN2011 software. 585 Cell area and cell number quantification 586 Cell area and cell number were determined by analyzing immunofluorescence images of 587 F-actin and nuclei via an custom-written image processing algorithm in MATLAB.

588
Binary masks of nuclei were generated by thresholding grayscale nucleus images; objects 589 in the binary mask were counted to determine total cell number. To determine cell size, 590 the centroid of each object in the binary mask was determined using the regionprops 591 function. Nuclei centroids were used to generate a Voronoi diagram, which consists of a 592 series of polygons that have edges that are equidistant from neighboring nuclei.

593
Previous studies have demonstrated that Voronoi diagrams reasonably predict cell 594 boundaries in an epithelial monolayer [40], and provide a more consistent quantification 595 of cellular size as opposed to quantification of protein markers in the cell-cell junction, 596 whose expression and localization changes as TGF-dose increases. Cell area was 597 calculated for each cell by summing the pixels in each Voronoi polygon, and were 598 averaged across the 250 µm x 250 µm colony. Spatial localization of cell number and 599 cell area were determined by binning nucleus centroids into a 5 x 5 grid. Cell counts in 600 each bin were totaled, and cell areas for each bin were averaged if the nuclei centroid 601 was contained within the bin. Spatial localization data was further combined into either 602 corner bins, edge bins, or interior bins, such that there were no overlap between the 603 three regions (i.e., corner bins were not included in the edge region).

605
To measure force on cell-cell junctions, Fluorescence Resonance Energy Transfer 606 (FRET)-based, full-length E-cadherin tension biosensors were stably transfected into 607 MDCK II cells. Epithelial square islands were cultured as stated above, and images 608 were acquired on a Zeiss LSM 710 laser scanning microscope using ZEN2011 software. 609 Briefly, mTFP (donor) and mEYFP (acceptor) fluorophores were imaged utilizing 610 spectral unmixing at 458 nm excitation. The acquired intensity images were manually 611 masked through ImageJ. Background subtraction and removal of saturated pixels was 612 then performed via an image processing algorithm in Python as previously 613 described [41]. FRET ratio was determined by obtaining the acceptor/donor ratio and 614 multiplying with a binary mask of the junctions. This allowed for inspection of FRET 615 pixels of interest within outlined cell-cell junctions.

616
Statistical analysis 617 Simulated and experimental data was exported to Prism 8 (GraphPad Software Inc) for 618 analysis. Statistical significance, indicated by a p-value less than 0.05, was determined 619 by one-way ANOVA across each TGF-1 dosage, ratio of interaction energies, and/or 620 spatial localization. about the centroid of the cluster (green dot). Junction forces (blue arrows) balance the 630 net force imbalance for a given cell. The linear system is constructed from the 631 mechanical equilibrium matrix and junction symmetry matrix (right). The mechanical 632 equilibrium matrix is constructed from the connectivity of each cell given by the 633 adjacency matrix and by applying the force balancing principle. The junction symmetry 634 matrix requires each junction force across a cell-cell adhesion to be equal and opposite. 635 S1 Video. Single cell without proliferation. Simulated cell organization for the 636 single cell FMA model as shown in Figure 1A. Movie corresponds to simulation of 1000 637 Monte Carlo Steps.