Optimal shapes and stresses of adherent cells on patterned substrates

We investigate a continuum mechanical model for an adherent cell on two dimensional adhesive micropatterned substrates. The cell is modeled as an isotropic and homogeneous elastic material subject to uniform internal contractile stresses. The build-up of tension from cortical actin bundles at the cell periphery is incorporated by introducing an energy cost for bending of the cell boundary, resulting to a resistance to changes in local curvature. Integrin-based adhesions are modeled as harmonic springs, that pin the cell to adhesive patches of a predefined geometry. Using Monte Carlo simulations and analytical techniques we investigate the competing effects of bulk contractility and cortical bending rigidity in regulating cell shapes on non-adherent regions. We show that the crossover from convex to concave cell edges is controlled by the interplay between contractile stresses and boundary bending rigidity. In particular, the cell boundary becomes concave beyond a critical value of the contractile stress that is proportional to the cortical bending rigidity. Furthermore, the intracellular stresses are found largely concentrated at the concave edge of the cell. The model can be used to generate a cell-shape phase diagram for each specific adhesion geometry.


I. INTRODUCTION
Living cells actively probe physical cues in their environment via receptor-ligand adhesion complexes that link the actomyosin cytoskeleton to the extracellular matrix (ECM) [1]. The cellular microenvironment, comprising of the ECM and of neighboring cells, imposes specific boundary conditions that can regulate physiological processes such as cell differentiation, division and motility, as well as cell architecture and polarity [2]. Myosin motors generate contractile stresses in the actin cytoskeleton that are transmitted to the substrate by focal adhesions. The traction stresses exerted by the cells on the substrate are thus very sensitive to the stiffness of the substrate as well as to the adhesion geometry. Cell morphology in turn is directly affected by traction stresses through the tension that builds up in the actomyosin stress fibers. It has been shown that the substrate stiffness plays a crucial role in regulating the cell spread area, the magnitude of traction forces and the cell morphology [3][4][5][6]. Much less explored is the role of adhesion geometry in regulating the spatial distribution of cellular stresses. Micropatterning has emerged as a powerful tool to investigate the interplay of mechanics and cytoskeletal architecture in controlling cell morphology by specific tuning of the geometry of the adhesion sites [7]. When plated on small micropatterns, cells are unable to grow, thus showing high apoptotic rate [8]. Large adhesive patches, in contrast, favor cell spreading and promote the assembly of contractile stress fibers along the cell's perimeter [9]. These peripheral stress fibers interconnect focal adhesions and yield concave arcs of constant curvature in the nonadherent portions of the cell boundaries. In addition, traction forces tend to localize in regions of high curvature at the boundary [10,11]. The model proposed here allows to separately study the roles of cell contractility and mechanical properties of peripheral cell fibers in controlling cell shape. Future comparison with experiments where both quantities can be perturbed using pharmacological interventions [9,12] may provide a quantitative understanding of the relative importance of boundary and bulk properties in determining steady state cell shapes.
Various successful theoretical models of single and multi-cell mechanics have been proposed over the past decade that address the role of ECM elasticity in regulating cell behavior [13]. Previous work has addressed the interplay between cell mechanics and geometry by either focusing solely on the elasticity of the cell boundary [14,15] or by considering only the bulk of the cell, described via continuum mechanics [16][17][18][19], by a cellular Potts model [20], or as a polymer network [21]. These models highlight the competing roles of cell contractility and substrate stiffness in regulating polymorphic cell shapes.
Continuum models of cell mechanics have assumed that the material constants describing the cell are spatially homogeneous. Cell material properties are, however, highly heterogeneous. In particular, experiments have shown strong differences in the mechanical properties of the bulk and boundary regions of the cell [22]. Increased tension and rigidity of cell boundaries can spontaneously arise during adhesion as a result of the assembly of peripheral stress fibers consisting of thin bundles of semiflexible actin filaments. Due to thermal and active forces these bundles considerably bend generating non-uniform peripheral tensions. Cell boundary can also resist changes in local curvature due to contact forces at the three-phase contact line between the cell, the substrate and the ambient medium. Motivated by these observations, in this paper we couple cell contour elasticity [14,15] to a continuum description of bulk cell mechanics [16,17] to investigate the cooperative roles of cortical elasticity, bulk elasticity and active contractility in controlling cell shapes on non-uniform adhesion pat-  terns. Non-uniform tension and elasticity is incorporated in the model by introducing a penalty for bending deformations of the cell periphery [15]. Using a combination of Monte Carlo simulation and analytical studies, we examine the interplay of bulk contractility and cortical tension in controlling morphological transitions in adherent cells and propose a cell shape phase diagram for specific adhesion geometries. An example of such a phase diagram for a cross-shaped adhesion pattern is shown in Fig. 1. See section III A for details.
The paper is organized as follows. In section II, we describe a continuum mechanical model for a thin adherent cell as an isotropic and homogeneous elastic material, subject to a homogeneous negative pressure, embodying active contractility. Cell-ECM adhesions are modeled as linear springs distributed non-uniformly along the cellsubstrate interface. Cortical tension is described via a penalty for bending deformations of the cell periphery. In Section III we discuss the steady shapes of cells adherent to different concave micropatterns obtained via Monte Carlo simulations (numerical details are given in the Appendix). Our simulations suggest a transition between convex and concave morphologies as a function of cell contractility and bending rigidity, that is captured by an analytically solvable model for the cell boundary presented in Section IV. We conclude with a brief discussion.

II. CONTINUUM MECHANICAL MODEL
We consider the mechanical equilibrium of a stationary cell strongly adherent to a soft elastic substrate. We assume that the cell's interior can be described as an isotropic homogeneous elastic material and neglect all dissipation. We further neglect all out-of-plane deformations of the cell and assume that its thickness is uniform throughout its entire area and remains unaffected by the substrate-induced deformations. The bulk elastic energy of our model cell is given by with E the three-dimensional Young's modulus, h the cell thickness, ν the Poisson's ratio, and u αβ = the strain tensor. We retain the nonlinearity in strain tensor to allow for the possibility of large strains that can arise even for small displacements [23]. The nonlinear terms essentially describe strain stiffening which is indeed expected in crosslinked actin networks [24]. The twodimensional displacement vector u is defined as u ( r 0 ) = r − r 0 , where r 0 ( r) is a material point before (after) the deformation. The integral is calculated over the area A 0 of the undeformed (reference) state and summation over pairs of repeated indices is assumed. Cell contractility arising from myosin motors is modeled as a homogeneous negative pressure resulting in an additional contribution to the cell's energy, given by where σ a > 0 is a parameter controlling active contractility, determine by concentration of myosin motors and rate of ATP consumption. At the continuum scale, it controls active contributions to the cell's surface tension [25]. Cell adhesion to the substrate is modeled via a harmonic potential with a position-dependent rigidity parameter Γ( r) with r a the position of focal adhesions on the substrate. The rigidity parameter Γ is nonzero over the adhesion region and zero elsewhere. Thus we allow for nonuniformity in the geometry of cell-substrate adhesions as can be realized experimentally using micropatterning techniques [9]. The assumption of local adhesive interactions with the underlying substrate strictly holds for elastic substrates that are much thinner than the cell perimeter or on soft microposts [26]. The rigidity parameter Γ depends on the elastic modulus of the underlying substrate as well as on the stiffness k f of focal adhesions.
For an elastic substrate of shear modulus µ s and thickness h s , with focal adhesion density ρ f , Γ is given by, therefore given by, . Finally, we assign a bending penalty to the cell's perimeter, reflecting the resistance of cortical actin bundles to changes in curvature, where κ is the bending rigidity, c = |γ (s)| is the curvature of the boundary, with γ (s) a parametric curve describing the cell boundary, and the line integral is calculated along the cell boundary. The optimal shape of the cell is obtained by minimizing the total mechanical energy E, that is given as the sum of elastic, active, adhesion, and boundary bending energies,

III. NUMERICAL SIMULATIONS
The minimal energy cell shapes have been determined numerically by a Monte Carlo study of a discrete representation of the continuum model introduced in Section II. The discrete representation of the undeformed cell is a triangulated disk. The initial configuration is built by randomly placing N ≈ 10 4 particles on a disk of radius R 0 . Particles are assumed to interact pairwise via a Weeks-Chandler-Andersen potential [27] and their positions are equilibrated using a standard Monte Carlo simulation with canonical (N V T ) Metropolis algorithm. A typical equilibrated configuration is stored and the positions of the particles of that configuration are then used as nodes to construct a Delaunay triangulation. The resulting triangulation for a V-pattern is shown in Fig. 2. We note that the initial density of points in the disk is chosen such that even in the equilibrium state there is a substantial overlap between neighbors, thus ensuring a densely packed distribution of points. As a result the equilibrium distribution of the interparticle distances is rather narrow and its mean, denoted as , represents a suitable unit of length. In the following, all distances are measured in units of and all energies are measured in units of Y 2 , where Y = Eh is the two-dimensional Young's modulus. The substrate rigidity Γ has units of Y 2 and the bending rigidity of cortical stress fibers has units of Y 3 . The low energy configurations are obtained using simulated annealing Monte Carlo (see Appendix for details).

A. Optimal shapes
We performed a series of simulations for V , U , and cross shaped micro-patterns, corresponding to the white dashed outlines shown in Fig. 3. The relaxed shapes for three non-convex patterns (U , V , and cross) obtained for a fixed value of the rigidity of adhesions Γ grey = 10 6 , and ν = 1/3 are shown in Fig. 3. In all cases, the non-adherent cell edges spanning two pinning regions are clearly concave. The relaxed shapes can be qualitatively compared with experiments on concave micropatterns [9]. Our simulations suggest that there is a transition between the concave and convex morphologies as a function of σ a and κ. In Figure 1, we show a sample cell shape phase diagram as a function of the cortical bending rigidity κ and the active contractile stress σ a for the cross shaped micropattern. The figure indicates that the cell boundary is concave at high values of the contractile stress σ a , whereas convexity is ensured at high values of bending rigidity. The phase boundary between convex and concave shapes appears to be linear in the σ a − κ plane. To justify this observation, in the next section we study analytically the shape of the cell boundary in the nonadhesive regions, considering small deformation about a circular configuration.

B. Optimal Stresses
Experiments probing the distribution of traction force density T ( r) exerted by cells adhering to soft substrates consistently show that such stresses are concentrated at the cell edges, and strongest in region of high cell curvature. Force balance requires T α = ∂ β σ αβ , where σ αβ is the two-dimensional stress tensor of the bulk cellular material, given by, The distribution of such internal stresses can therefore be inferred experimentally from traction force microscopy measurements [28]. Internal stresses of adhering cells are found to be concentrated at the cell's interior, with a maximum value proportional to the active cell contractility, here σ a . To highlight the role of patterned adhesion on the spatial distribution of cellular stress, we display in Fig. 4 the spatial distribution of the so-called Lamé's stress ellipses [29] for the elastic part σ el αβ of the stress tensor. The constant active contribution σ a has been subtracted out to highlight spatial variations. As a result, the displayed stress is largest at the cell edges. The Lamé's stress ellipses are obtained by computing the elastic part of the two-dimensional stress tensor at a representative subset of triangles in the Delaunay grid. This is achieved by directly evaluating the expression in eqn (5), excluding the active term. We then compute the low and the high eigenvalues, σ min and σ max , respectively, of the elastic stress tensor of a given triangle. Note that since the stress tensor is symmetric, its eigenvalues are always real. The length of the major and minor semiaxes of each ellipse are then given by σ max and σ min , respectively, whereas the orientation of the ellipse axes is given by the directions of the corresponding eigenvectors. As expected, and consistent with experiments [7], elastic stresses are concentrated at the free boundaries of the adherent cell. Boundary stresses along free edges connecting two adhesion points are directed normal to the edge whereas they are oriented along the edge near the adhesion points. This is most evident in the cross-shaped pattern, Fig. 4(c). The large stresses in the convex regions of the cell spilling outside straight portions of the pinning regions (see, for instance, the V-shape pattern, Fig. 4(a)) are largely an artifact of our model. They arise because we have introduced excluded volume interactions to prevent self-intersections of the triangulation. In other words, we assign a hard core radius of 0.25 to each vertex of the triangulation, such that no two vertices can came closer than 0.5 from each other. Once this limit has been reached, the excluded volume prevents further collapse of the cell, thus accounting for the presence of a sizeable portion of the cell that extends outside the pinning region. While steric effects are present in vivo and may describe for instance the role of structural elements capable of carrying compressive loads, such as microtubules, cells on synthetic substrates generally almost completely conform to the micropattern by changing their thickness. This is not possible in our strictly two-dimensional model. As a result, the model captures well the behavior of "free" cell edges spanning two adhesion points, but has limitations for describing the behavior of cell boundaries along straight pinning regions.

A. Strong pinning at adhesions
For small deformations about an initially circular configuration of radius R 0 , the cell boundary in the nonadhesive region can be parametrized using polar coordinates as, r(θ) = R 0 + u r (θ), where u r is the radial component of the displacement field at the cell boundary, u = (u r , u θ ). Thus, u is solely a function of the angular coordinate θ. In mechanical equilibrium, boundary force balance along the normal and tangential directions requires where s is the arc-length parameter and n and t are unit vectors normal and tangent, respectively, to the unperturbed cell boundary. Tangential force-balance in polar coordinates reduces to σ rθ = 0, which leads to the relation u θ = ∂ θ u r . Thus, the normal component of the elastic stress is given by σ el rr λ∂ θ u θ /R 0 = λ∂ 2 θ u r /R 0 , where λ = Y ν/(1 − ν 2 ) is the Lamé elastic constant. 0 /κ = 1 (solid, black), σaR 3 0 /κ = 10 (dashed, blue) and σaR 3 0 /κ = 40 (dotted, red), where 0 < θ < π/2, corresponding to a cross pattern of width zero. (c) Shape phase diagram for the contractile string pinned to a cross micropattern of width w, obtained from the solution to eqn (8). Bending rigidity κ and contractile stress σa are given in units of Y 3 and Y respectively, corresponding to the parameters in Fig. 1, where R0/ = 50, w/ = 20, and ν = 1/3.
Furthermore, for small deformations u r , the boundary curvature can be expanded as, Using eqns (6a)- (7), and letting ds = R 0 dθ, we obtain an equation for the boundary profile, whereũ = u r /R 0 . Without loss of generality, we can consider solution in the interval θ ∈ [0, φ], where φ is the angular width of the nonadherent region and depends on the geometry of the adhesion pattern. The boundary conditions for the case of strong pinning at adhesions are given by:ũ(0) =ũ(φ) = ∂ θũ (0) = ∂ θũ (φ) = 0. The full solution of eqn (8) is analytically tractable but cumbersome. We instead discuss the solutions in two limiting cases in terms of the dimensionless parameter K = 2κ/λR 3 0 , reflecting the relative contributions of bending and bulk elasticity. This parameter can also be written as K = (ξ/R 0 ) 3 in terms of the ratio of a length scale ξ = [2κ/λ] 1/3 to the undeformed cell radius R 0 . The length scale ξ described the interplay between bulk elasticity and boundary tension in controlling the response of the cell. When ξ R 0 (corresponding to K 1) the cel deformation is controlled by the cortical tension at the boundary and the curvature is given by Bulk elasticity drops out and the behavior is controlled by the ratio σ a R 3 0 /κ of contractility to bending rigidity. The curvature has a minimum at the center of the nonadherent segment. Thus, as one increases contrac-tility σ a , a region of negative curvature develops near θ = φ/2, which grows upon increasing σ a until convexity is retained within a small neighborhood of the adhesion patch. The onset of concavity is thus given by the condition of reality to the solution of c(θ) = 0, which gives the condition, Since concave shapes are commonly observed in experiments [9], we now turn to estimate the critical value of σ a as predicted by our model in order to compare it with experimentally reported values for σ a . The bending rigidity of cortical stress fibers can be estimated as κ ∼ π 4 E act r 4 s , where E act is the Young's modulus of actin and r s is the typical radius of the stress fibers. Using E act 2.6 GPa [30] and r s ∼ 0.1 µm, we get κ ∼ 2.0 × 10 −19 Nm 2 . Using this value in eqn (10) for φ ∼ 2π/3, corresponding to a thin V-pattern, we get value for the critical σ a ∼2.7 nN/µm. This is indeed the order of magnitude value for active stress or surface tension reported in experiments for adherent epithelial cells on continuous elastic substrates or endothelial cells on microposts [14,25]. In the opposite limit of K 1 the deformation is controlled by bulk elasticity and the curvature is given by The condition of concavity is given by σ a > λ 1 φ 2 /8−1 . In the case when ξ is comparable to R 0 , a simpler solution of eqn (8) can be obtained by neglecting the fourth order gradient term and also the derivative boundary conditions. The crossover to concave profiles can be approximated by the following interpolating form between the two limiting cases, In the general case of eqn (8), the solution for curvature and the radial displacement is given in Fig. 5a,b for three different values of σ a R 3 0 /κ that compares the relative strengths of contractility to bending deformations. Furthermore, to compare numerically with the simulation results for the shape phase diagram of the adherent cell, we show the concave-convex phase boundary in σ a − κ plane in Fig. 5c, for a cross-shaped micropattern using the same parameters as used in Fig. 1. The resultant phase diagram is in good order-of-magnitude agreement with the simulation results, and the discrepancy in numerical values possibly arise from neglecting non-local bulk elasticity in the theoretical analysis.

B. Soft pinning
We now consider the case of soft pinning, where the free cell boundary is anchored to soft springs at the adhesion sites. Equation (8) is now solved with the boundary conditionsũ(0) =ũ(φ) = δ and ∂ θũ (0) = ∂ θũ (φ) = 0, where we have introduced an unknown displacement δ of the ends of the segment, which can be self consistently determined by minimizing the total energy of the deformed configuration with respect to δ. For simplicity we ignore bulk elasticity and consider the limit K 1. The total energy of the deformed configuration is then given by U = κR 0 φ 0 dθc(θ) 2 + k s δ 2 , where k s = ΓA f and A f is the cross-sectional area of focal adhesions. Note that the contribution due to contractility vanishes in the final energy due to the derivative boundary conditions oñ u. The onset of concavity now depends on the substrate stiffness k s and the condition for convex-concave transition is given by (13) Thus, stiffer adhesions with k s R 3 0 κ, promote concavity transition at a much higher value of contractility. It is favorable for a cell to invaginate at the free edges if the anchoring at adhesions is softer than the effective bending stiffness κ/R 3 0 .

V. CONCLUDING REMARKS
Using a simple continuum model, coupling bulk and contour mechanics, we investigate the equilibrium shapes and stresses of adherent cells on substrates with various adhesion patterns. A continuum model without contour elasticity have been studied previously by two of us on convex patterns [16], which was successful in capturing distribution of traction and cellular stresses and their dependence on substrate physical properties [26]. Here we focus on the shape and geometry induced stresses of non-adherent cell edges on concave micropatterns. We demonstrate numerically and analytically that the curvature of the non-adherent cell boundary can undergo a shape transition from convex to concave morphology, controlled by the interplay of contractility and bending rigidity. Stiff boundaries with low contractility relax to convex shapes, whereas at higher values of contractility, non-adherent cell edges attain a concave morphology. Previous work has shown that contractile cable network models are capable of reproducing the invaginated circular arc morphology of cell edges connecting strongly adhering sites [21]. Here we demonstrate that simple continuum whole-cell models can also predict qualitatively cell shape and the transition between convex and concave cell edges, provided a bending rigidity describing cortical tension is included. For parameters realistic to experiments (see section 4.1) our model suggest that cells prefer to invaginate at their free edges, such that the effective boundary stiffness on non-adhesive zones are softer than myosin induced contractile stresses. Images of actin from experiments on concave micropatterns do indeed show the formation of long and thin stress fibers that are invaginated on non-adherent edges [9], indicating a softer cortical rigidity. In addition, elastic stresses are found to be higher along the free cell boundaries than in the neighborhood of adhesions, since in the absence of mechanotransduction cellular forces along free edges are not shared by the substrate. Previous theoretical study with only contour elasticity indicated that substrate stiffness and contractility can cooperatively control cell morphology and induce hysteresis at the onset of convex-concave transition [15]. Here we show that even in the presence of rigid adhesions, cell shape can be controlled by regulating the cortical bending rigidity and contractility. Bending rigidity can be experimentally controlled by regulating the amount of actin cross-linking proteins that can impact stress fiber thickness and rigidity, whereas myosin based contractility can be perturbed using the conventional inhibitor Blebbistatin. One limitation of our model is that it is strictly two-dimensional and does not allow for changes in the cell thickness. Due to the presence of steric interactions in the finite element simulations, the cell edges on flat adhesive segments do not fully relax to the flat morphology, but maintain a convex shape. This is in contrast to real cells that contract to adjust to the shape of the micropattern. A fully three-dimensional model can overcome this difficulty, and is a natural extension of our present work.

A. Discrete Model
The discrete version of the elastic and active contraction energies can be expressed as a sum over triangles of the triangulation,