Introduction

Small self-assembled monolayers are common in chemistry, condensed matter physics, and biology.1,2,3 These monolayers often exhibit rich morphologies, arising from a combination of short-range inter-monomer interactions and long-range interactions through elastic or electrical fields. For example, a graphene or MoS2 monolayer on a substrate may form wrinkles with different wavelengths due to the elastic, edge, and interfacial effects.4,5 These effects may also present in multicellular epithelial ensembles. Yet, different from atomic monolayers or thin films in condensed matters, a multicellular monolayer is biochemically active in that it probes and senses physical cues from its surroundings and responds by a cascade of biochemical signals.6,7,8 The biochemical activities are accompanied by mechanical force generation and transmission,9,10 diffusion and assembling of mechanical force responsive molecules,11,12 and morphogenesis of the entire multicellular monolayer.13,14 Here we address the cellular force landscapes of equilibrated multicellular epithelial monolayers and its implications to cell fate and functions.

We have examined the HCT-8 epithelial cell colonies (a colony consists of hundreds of cells packed together through cell–cell adhesion.), as an example of cohesive multicellular monolayers. Upon seeding adherent HCT-8 cells onto polyacrylamide (PAA) gels in vitro, they proliferate and aggregate to form cohesive multicellular monolayers of different sizes and morphologies (Fig. S1). The proliferation of the cells slows down and the colonies size and shape nearly unchanged beyond 3 days of culture. In the course of morphogenesis, mechanosensation of mechanical cues from the substrate through cell surface integrin receptors activates a complex biochemical-signaling network, leading to increased actomyosin contractility.15 The increased intracellular contraction feeds back to the mechanotransduction pathway by promoting integrin clustering,16,17,18 focal adhesion assembling,12,19 and cell–cell adhesion and cytoskeletal remodeling.20 As cell–cell adhesion junctions and focal adhesions are both biochemically and structurally linked, intracellular contraction invokes cross-talks between them,21,22 generating extracellular traction onto the substrate and intercellular tension in a coordinated fashion that are counter-balanced by cell body stress. We have therefore developed an energy functional, combining monolayer and substrate elasticity, integrin-mediated interfacial adhesion, and chemical potential of adhesion molecules, to determine the cellular force landscapes of the cohesive multicellular monolayers.

Results

To a minimal setting, we consider an epithelial monolayer with a total of N cells seeded on PAA gels coated with fibronectin to which integrins bind and form focal adhesion points (see Fig. 1). The extracellular traction and intercellular tension originated from actomyosin contraction displace both the substrate and the monolayer. Adopting a continuum view within the framework of elasticity, we denote the displacement fields of the monolayer and the substrate by u and \(\widehat {\bf{u}}\), and the respective conjugate stress fields by σ and \(\hat \sigma\). Focal adhesion points furnish the structural units that transmit forces between the monolayer and the substrate. The bridging role and spatially heterogeneous distribution of focal adhesion points warrant a microstructure-based description. To capture the molecular basis of focal adhesion formation through integrin clustering, we categorize the integrin receptors on the cell membrane into two distinct phases: a freely diffusive phase with a density of ϕi,f (number per unit area) and a bound phase to the ligands on the substrate with a density of ϕi,b for the i-th cell. Conservation of the receptors expressed on the membrane of each cell requires \({\int} {_{\Omega _i}} (\phi _{i,b} + \phi _{i,f}){\mathrm{d}}\Omega = \phi _{i,0}\Omega _i\) with ϕi,0 being a constant for the i-th cell that covers an area of Ω i . These two phases, along with a vacancy phase of density ϕi,v, occupy all the molecular sites on the cell membrane. The three phases thus conserve at any sites on the cell membrane: ϕi,f + ϕi,b + ϕ i,v  = 1. Only the bound phase constitutes the focal adhesion points and transmits the force. Within the focal adhesion points, each receptor-ligand pair is modeled as a linear spring with a spring constant k. The stretching force in the pair is F = kΔu, where \(\Delta {\bf{u}} \equiv \widehat {\bf{u}} - {\bf{u}}\). Accordingly, the traction force generated by the receptor-ligand pair and transmitted to the substrate is T = n l ϕ b F, where n l is the lattice site per unit area on the cell membrane and ϕ b  = {ϕi,b} lumps the spatially organized focal adhesion points of all the constituent cells.

Fig. 1
figure 1

Schematics of force generation mediated by focal adhesion formation through integrin clustering in a multicellular monolayer. a Randomly distributed, freely diffusible integrin (red dots) on cell membrane. b Focal adhesion formation through integrin clustering upon cell-substrate interactions. c Molecular structure of the focal adhesions

Full size image

With the above settings, the free energy functional of the cohesive monolayer takes the following form:

$$\begin{array}{*{20}{l}} {{\mathrm{II}}\left( {\phi _{{\it{i}}{\mathrm{,}}{\it{b}}}{\mathrm{;}}{\bf{u}}{\mathrm{;}}\widehat {\bf{u}}} \right)} \hfill & = \hfill & {\mathop {\sum}\nolimits_{i=1}^N {W_i\left( {\phi _{i,f};\phi _{i,b};\phi _{i,v};\Delta {\bf{u}}} \right)} } \hfill \\ {} \hfill & {} \hfill & { + \frac{h}{2}{\int}_\Omega {\sigma :\varepsilon {\mathrm{d}}\Omega } } \hfill \\ {} \hfill & {} \hfill & { + {\int}_{\partial \Omega } {{\mathrm{\Gamma d}}\ell } } \hfill \\ {} \hfill & {} \hfill & { + \frac{1}{2}{\int}_V {\hat \sigma } :\hat \varepsilon {\mathrm{d}}{\it{V}},} \hfill \end{array}$$
(1)

where the first term is the sum of the free energies associated with the receptors of all the constituent cells, and the next three terms are constructed on the continuum level without differentiating individual cells. This treatment features a multiscale hybrid form of the thermodynamic model. In Eq. (1), the first term denotes the chemical energy, which includes the chemical energy in the reference configuration, the mixing entropy of the three phases, the gradient energy of three different phases, the stretch energy stored in each ligand-receptor bond, and chemical energy release upon the binding of a receptor-ligand pair, as detailed in Supplementary Information. The second term is the strain energy of the cell sheet, where h and Ω are the thickness and area of the monolayer. The third term is the line tension energy on the boundary of the cell sheet, where Γ is the line tension energy density (per unit length). The last term is the strain energy stored in the substrate. Minimizing the free energy functional with respect to its independent variables yields the partition of the receptors in each cell, as:

$${\mathrm{ln}}\left( {\phi _{{\it{i}}{\mathrm{,}}{\it{b}}}{\mathrm{/}}\phi _{i,f}} \right) - \beta \nabla ^2\left( {\phi _{i,b} - \phi _{i,f}} \right) = \gamma {\mathrm{/}}2 - \Delta \mu ^0,i = 1,...,N,$$
(2)

where β is the coefficient of the gradient energy, Δμ0 is the difference of the reference chemical potentials of the bound and free receptors, and γ is the chemical energy release upon the binding of a receptor-ligand pair. The chemical energy release is equal to the work done by the force that stretches the receptor-ligand bond, γ = F Δu. One notes from Eq. (2) that ϕi,b increases with tension, manifesting that the focal adhesions are mechanical responsive and tension promotes focal adhesion formation.12,23

Minimizing the free energy functional with respect to the displacement field u gives rise to the mechanical equilibrium condition of the cell sheet

$$\nabla \cdot \sigma + Y\Delta {\bf{u}}/{\it{h}} = 0\,{\mathrm{in}}\,\Omega$$
(3)

where Y = n l kϕ b characterizes the coupling between the monolayer and the substrate. The term YΔu is the force sustained in the focal adhesion points and transmitted to both the monolayer and the substrate. Owing to the thinness of the monolayer, this force is treated as a body force that is uniformly distributed over the monolayer thickness. The monolayer is subjected to the traction boundary condition at its periphery, as

$$\sigma \cdot {\bf{n}} = - \frac{{\Gamma \kappa}}{h}{\bf{n}}$$
(4)

where n and κ are the outer unit normal vector and the local curvature of the monolayer boundary, respectively. We model the epithelial monolayer as a two-dimensional (2D) isotropic, elastic sheet, with Young’s modulus E and Poisson’s ratio ν. The intracellular contraction is modeled by an active pressure σ A within the epithelium, resembling a force induced by thermal cooling. The constitutive relation of the monolayer is set as,24

$$\sigma _{ij} = \frac{E}{{1 - v^2}}\left[ {v\delta _{ij}\nabla {\bf{u}} + \frac{{1 - v}}{2}(u_{i,j} + u_{j,i})} \right] + \sigma _A\delta _{ij}$$
(5)

Similarly, the substrate satisfies the mechanical equilibrium \((\nabla \cdot \hat \sigma = 0)\) and surface traction boundary condition (T = |YΔu). The substrate is modeled as an elastic half space with a Young’s modulus \(\hat E\) and Poisson’s ratio \(\hat v\).

We employ finite element methods to solve the coupled governing equations. We use cubic triangular elements to discretize the geometry of the cell colonies. The meshes used here are fine enough to follow the curvature of the boundaries of the cell colonies and to accurately interpolate the displacement and stress inside the colonies.

Active, nonlinear cell-matrix coupling

To elucidate the active coupling between the multicellular monolayer and the substrate, we note that the intracellular contractility generates in the monolayer such a displacement field that monotonically decreases in magnitude (|Δu|) from the periphery to the center (Fig. S1), analogous to the in-plane shrinking of a thin film under thermal cooling. The intracellular contraction is transmitted through the ligand-receptor pairs, generating traction force onto the substrate. The traction force spatially varies with the displacement field (|T| ~ | Δu), and is expected to be high at the periphery of the monolayer because of large |Δu| therein. Ligand-receptor pairs at the periphery of the monolayer are also under high tension since |F| ~ | Δu|, which promotes integrin clustering (increasing ϕ b ) and stabilizes focal adhesion points,12 as indicated by Eq. (2). The large ϕ b further strengthens traction force as |T| ~ ϕ b . The traction force field then feeds back to the monolayer, establishing the intercellular tension and cell body stress. Thus, the thermodynamic model faithfully follows the positive feedback mechanical network in force generation and transmission. As shown in Fig. 2a,b, the active monolayer-substrate interaction results in co-localization of focal adhesions and traction force of the HCT-8 monolayers, obtained by the immunofluorescence of the focal adhesion kinase (FAK) (Fig. S3) and traction force microscopy (TFM) (Fig. S4), respectively. The traction force and focal adhesions are generally localized at the boundary layer of ~50 μm, a width that spans roughly 5 HCT-8 cells.

Fig. 2
figure 2

Nonlinear reciprocal cell-matrix coupling. Left: a Immunofluorescence staining shows localization of focal adhesion kinase (FAK) at the boundary layer of a monolayer. b The traction force map measured by TFM. c, d The traction force maps predicted by the previous uniform coupling model (c, constant Y)24,25 and the current nonlinear coupling model (d, nonlinear Y) model, subtracted by the map from TFM in (b). The darker the pixel, the better the agreement of the traction force. e The decay of the radial component of traction force from the periphery to the centroid of HCT-8 colony of 186 μm in radius, shown in (b), predicted by the models and TFM. All the results are normalized by the maximum average radial traction forces

Full size image

The coupling strength between the monolayer and the substrate, characterized by Y, is highly nonlinear (Fig. S2). The nonlinearity arises from the mechanically responsiveness of the focal adhesion points, manifested by the nonlinear dependence of ϕ b on the stretch |Δu|. The coupling strength is thus spatially heterogeneous, in contrast to the assumption of uniform coupling (constant Y in space) in previous models.24,25 The characteristics of monolayer-substrate coupling dictate traction force landscapes, as theoretically examined below by considering a circular shaped monolayer with radius R. For uniform coupling, the radial displacement of the circular monolayer decays with a modified Bessel function, u r  ~ I1(r/R), so as the traction force.24 In contrast, for the nonlinear coupling derived here, the radial displacement in a first-order approximation decays with of the derivative Airy function of the second kind, u r  ~ (R/r)A 2 '(r/R), and the traction force with A 2 '(r/R) (see SI text). This decay is markedly faster than that under uniform coupling, suggesting a stronger localization of the traction force and focal adhesion points at the boundary layer. Figure 2c, d depict the traction force maps of an irregularly shape monolayer, predicted by the uniform coupling model (c) and the current nonlinear model (d). For a better comparison with TFM, the maps are color-subtracted by that obtained from TFM in Fig. 2b. A darker pixel color means a better agreement with the TFM data. The difference of the models is more clearly seen in Fig. 2e. Here each data point, either from TFM or model predictions, is taken as the average over the small region (r − δ,r + δ), where r measures from the centroid of the monolayer, and δ << R. The nonlinear model (solid black line) yields a faster decay of the traction force from the periphery to the centroid of the monolayer than the uniform coupling model (dashed green line), which is consistent with the theoretical analysis for circular shaped monolayers and better compares to the experimental results from TFM (open circles).

Intracellular contractility depends on substrate stiffness

A contractile cell is able to sense the substrate stiffness and responds by modulating its contractility and cytoskeletal remodeling.26 The active stress σ A , a measure of the contractility, is thus not a constant but depends on the substrate stiffness. To determine the functional dependence of the active stress on the substrate stiffness, we examine the HCT-8 monolayers of roughly the same radius R ~ 150 μm but with hydrogels of different stiffness. Here the equivalent radius of the monolayer is given by \(R = \sqrt {A/\pi }\), and A is the area of the monolayer. Gels with stiffness larger than 50 kPa are physiologically irrelevant7 and thus not considered here. The traction force profiles of the monolayers are then measured by TFM. The traction forces are also predicted by the thermodynamic model with assumed active stress. The active stress for different gel stiffness is then determined by matching the numerically obtained traction force profile over the boundary layer with the experimental measurements. Our data matching shows that the active stress monotonically increases with the gel stiffness and saturates when the substrate stiffness reaches a critical value27,28 (\(\hat E > 30{\mathrm{kPa}}\), Fig. 3). The saturation in the active stress manifests the limit of cell internal machinery in driving crosslinking actomyosin motors.

Fig. 3
figure 3

Substrate stiffness mediates cell contractility. Active stress vs. Young’s modulus of substrate \(\hat E\) obtained by experimental-numerical matching of the traction forces

Full size image

Intracellular traction and intercellular tension are both gel stiffness and colony size dependent

With the set of governing equations and physiologically relevant parameters (see Supplementary Information), we set forth to simulate the traction force profiles of the monolayers with different sizes and on substrates of different stiffness (with parameters in the supplementary material), as shown in Fig. 4a1, a2 for the phase contrast images. The traction force map of the colonies is first measured by TFM (Fig. 4b1, b2). The morphologies of monolayers are digitalized and used as the reference configurations for modeling, with special attentions paid to the details at the periphery of the monolayers. Our model captures the overall traction map remarkably well, as shown in Fig. 4c1, c2. Traction forces are highly localized at the boundary layer of the monolayers, consistent with the localized distribution of the focal adhesion points. The model predicts that the traction force is gel stiffness dependent (Fig. 4c1), originated from the stiffness-dependent cell contractility. The predicted average traction forces are 61.7 kPa, 105.5 kPa, and 138.6 Pa for the colonies on 4.5 kPa, 20.7 kPa, and 47.1 kPa gels, respectively. Here the traction forces are averaged over the boundary layers of ~50 μm in width. At the periphery of the colonies, the traction force is markedly larger at convex regions than concave regions, owing to the curvature-dependent surface tension. This follows that the traction force is also colony size dependent since the average curvature decreases with the colony size (Fig. 4c2). Specifically, the average traction forces within the 50 μm thick boundary layer are 122 Pa, 110 Pa, and 105.5 Pa for 54 μm, 107 μm 186 μm colonies, respectively. The model predictions on the monolayer size and gel stiffness effects agree very well with the TFM.

Fig. 4
figure 4

Traction force is both substrate stiffness (Left panel) and colony size (Right panel) dependent. a1, a2 Phase contrast images of HCT-8 cell colonies on different gel stiffness (a1,with gel stiffness of 4.5 kPa, 20.7 kPa and 47.1 kPa, respectively; scale bar: 90 μm) and with different monolayer sizes (a2, scale bar: 45 μm) on 20.7 kPa gels. b, c The traction force profiles of the cell colonies obtained by TFM (b1, b2) and modeling (c1, c2), respectively. d,e The intercellular tension maps of the cell colonies obtained by MSM (d1, d2) and modeling (e1, e2)

Full size image

Discussion

Cell contractility is partially transmitted to the neighboring cells, generating intercellular tension and the cell body stress. Pertaining to the 2D thin film structure of the monolayer, we here plot the 2D hydrostatic tension (the average of two principal stress components) in the monolayer. The cellular stress can be obtained from the monolayer stress microscopy (MSM, see SI text)29,30 as well, where the traction force map obtained from TFM is used as a boundary condition for finite element analysis (Fig. 4d1, d2). As a general trend, the cellular stress is minimal at the periphery of the cell sheet, ramps up over the boundary layer, and finally reaches a uniform tension. The cellular stress is spatially oppositely correlated with the traction force. This opposite correlation stems from the crosstalk between the focal adhesions and intercellular adhesion junctions in transmitting the intracellular contraction. At the boundary layer, cell contraction is mostly transmitted through the focal adhesions to generate extracellular traction. Whereas at the interior of the cell sheet, cell contraction is mostly transmitted to the neighboring cells through intercellular adherent junctions, generating uniform intercellular tension and cell stress. Our model predicts that cellular stress scales with the gel stiffness (Fig. 4e1) and the colony size (Fig. 4e2) in the similar manner as the traction force at the boundary layer of the colonies. That is, the stiffer the gels and the smaller the colonies, the larger the 2D hydrostatic tension at the interior of the colonies, consistent with MSM results (Fig. 4d1, d2).

In conclusion, we have developed a molecularly based thermodynamic model that simultaneously incorporates the elastic, edge, and interfacial effects in cohesive epithelial monolayers. Consistent to the mechanosensation of focal adhesions and subsequent activation of the biochemical signaling cascade, the model captures the active, nonlinear cell-matrix coupling. With these features, the model predicts monolayer size and substrate stiffness dependent traction force and intercellular tension, validated by microscopy and immunofluorescence studies. These results suggest that substrate stiffness as a mechanical cue defines how efficient cellular tension and cell body stress can be generated, posing a controlling factor for tensional homeostasis.31 For extremely compliant or stiff substrates, abnormally low or high intercellular tension causes the disruption of stress homeostasis, resulting in apoptosis or phenotypic changes. Our results thus infer how cell fate can be steered by the local mechanical environment.

Our model has focused on the thermodynamics of cohesive epithelial monolayers constituted of slow-migrating epithelial cells. The active stress, treated here as a uniform pressure, can also be extended to be anisotropic to capture polarization and re-orientation of the cells under external chemical, mechanical and topographic stimuli. Further, our model sets a firm ground toward a versatile computational framework to address the kinetics of morphogenesis and collective migrations of monolayers, a more complex problem as the multicellular ensemble dynamically evolve due to cell proliferation and migration. Such a computational framework offers a cost-effective alternative to TFM for understanding the morphogenesis and disease of multicellular epithelia.

Methods

The weak form of the governing equations was developed and the finite element analysis was performed using the COMSOL Multiphysics package. Cubic triangular elements were used to interpolate the displacement and stress inside the colonies. The relevant parameters are listed in the Supplementary Information.

For the experimental methods and materials, PAA hydrogels of various stiffnesses were prepared by tuning the concentrations of acrylamide and bis-acrylamide. The PAA gel surface was coated with fibronectin for cell seeding. Fluorescent beads were embedded in a single plane beneath the hydrogel surface to track the displacement caused by cell traction. The HCT-8 cells were seeded at a density of 2000 cells/cm2 onto gels and incubated for 24–72 h, which allows cells to adhere, spread and grow into isolated multicellular colonies of different sizes. Three-dimensional (3D) images of cell colonies were taken by using a laser scanning confocal fluorescence microscope (Olympus FV10i, Japan) with a 60 × water-immersion lens (NA = 1.2, Olympus, Japan) to estimate the cell thickness. TFM and MSM were then applied to quantify traction force and intercellular force maps. Details can be found in the Supplementary Information.

Data availability statement

The authors declare that all the data supporting the findings of this study are available within the paper and its supplementary information files.