Mechanical interplay between cell shape and actin cytoskeleton organization

TL;DR

This study explores how the organization of the actin cytoskeleton influences cell shape and vice versa, using a combined analytical, simulation, and experimental approach to reveal a feedback mechanism governed by a dimensionless anchoring number.

Contribution

The paper introduces a novel model linking actin cytoskeleton organization with cell shape, incorporating a feedback mechanism controlled by the anchoring number, validated by experiments.

Findings

Cell edge shape approximated by elliptical arcs in aligned cytoskeletons

Eccentricity of cell shape reflects internal stress anisotropy

Model accurately predicts cell shape and cytoskeleton structure

Abstract

We investigate the mechanical interplay between the spatial organization of the actin cytoskeleton and the shape of animal cells adhering on micropillar arrays. Using a combination of analytical work, computer simulations and in vitro experiments, we demonstrate that the orientation of the stress fibers strongly influences the geometry of the cell edge. In the presence of a uniformly aligned cytoskeleton, the cell edge can be well approximated by elliptical arcs, whose eccentricity reflects the degree of anisotropy of the cell's internal stresses. Upon modeling the actin cytoskeleton as a nematic liquid crystal, we further show that the geometry of the cell edge feeds back on the organization of the stress fibers by altering the length scale at which these are confined. This feedback mechanism is controlled by a dimensionless number, the anchoring number, representing the relative weight of surface-anchoring and bulk-aligning torques. Our model allows to predict both cellular shape and the internal structure of the actin cytoskeleton and is in good quantitative agreement with experiments on fibroblastoid (GD$\beta$1,GD$\beta$3) and epithelioid (GE$\beta$1, GE$\beta$3) cells.

Equations

$${\xi }_t{\partial }_t\mathbf{r}={\partial }_s{\mathbf{F}}_{\mathrm{cortex}}+({\hat{\mathbf{\Sigma }}}_{\mathrm{out}}-{\hat{\mathbf{\Sigma }}}_{\mathrm{in}})\cdot \mathbf{N},$$

$${\hat{\mathbf{\Sigma }}}_{\mathrm{out}}-{\hat{\mathbf{\Sigma }}}_{\mathrm{in}}=\sigma \hat{\mathbf{I}}+\alpha \mathbf{n}\mathbf{n},$$

$${\xi }_t{\partial }_t\mathbf{r}={\partial }_s(\lambda \mathbf{T})+\sigma \mathbf{N}+\alpha (\mathbf{n}\cdot \mathbf{N})\mathbf{n}.$$

$$0=({\partial }_s\lambda )\mathbf{T}+(\lambda \kappa +\sigma )\mathbf{N}+\alpha (\mathbf{n}\cdot \mathbf{N})\mathbf{n},$$

$$({\partial }_s\lambda )\mathbf{T}+(\lambda \kappa +\sigma )\mathbf{N}=0.$$

$$\lambda \mathbf{T}+(\sigma \hat{\mathbf{I}}+\alpha \mathbf{n}\mathbf{n})\cdot {\mathbf{r}}^{\perp }={\mathbf{C}}_1,$$

$$\lambda \cos \theta$$

$$\lambda \sin \theta$$

$$\frac{1}{\gamma }(x-x_c{)}^2+(y-y_c{)}^2=C_2,$$

$$x_c=\frac{C_{1y}}{\sigma +\alpha },y_c=-\frac{C_{1x}}{\sigma },\gamma =\frac{\sigma }{\sigma +\alpha }.$$

$${\lambda }^2={\sigma }^2(y-y_c{)}^2+(\sigma +\alpha {)}^2(x-x_c{)}^2.$$

$$\frac{{\lambda }^2}{{\sigma }^2}=C_2\frac{1+{\tan }^2\theta }{1+\gamma {\tan }^2\theta }.$$

$$\lambda (\theta )={\lambda }_{\min }\sqrt{\frac{1+{\tan }^2\theta }{1+\gamma {\tan }^2\theta }}.$$

$$\frac{{\sigma }^2}{\gamma {\lambda }_{\min }^2}(x-x_c{)}^2+\frac{{\sigma }^2}{{\lambda }_{\min }^2}(y-y_c{)}^2=1.$$

$$\begin{gathered} \frac{{\sigma }^2}{{\lambda }_{\min }^2}{\left[(x-x_c)\cos {\theta }_{\mathrm{SF}}+(y-y_c)\sin {\theta }_{\mathrm{SF}}\right]}^2 \\ +\frac{{\sigma }^2}{\gamma {\lambda }_{\min }^2}{\left[-(x-x_c)\sin {\theta }_{\mathrm{SF}}+(y-y_c)\cos {\theta }_{\mathrm{SF}}\right]}^2=1. \end{gathered}$$

$$\kappa =-\frac{1}{\gamma b}{\left(\frac{1+\gamma {\tan }^2\theta }{1+{\tan }^2\theta }\right)}^{\frac{3}{2}},$$

$${\kappa }_{\min }$$

$${\kappa }_{\max }$$

$$\bm{F}_{0}=\lambda_{\mathrm{min}}\Biggr{[}\Biggr{(}\frac{d}{2b}\sin\phi+\frac{\rho}{b}\cos\phi\Biggr{)}\>\bm{n}^{\perp}\\ +\Biggr{(}-\frac{1}{\gamma}\frac{d}{2b}\cos\phi+\frac{\rho}{b}\sin\phi\Biggr{)}\;\bm{n}\,\Biggr{]},$$

$$\bm{F}_{1}=\lambda_{\mathrm{min}}\Biggr{[}\Biggr{(}\frac{d}{2b}\sin\phi-\frac{\rho}{b}\cos\phi\Biggr{)}\>\bm{n}^{\perp}\\ +\Biggr{(}-\frac{1}{\gamma}\frac{d}{2b}\cos\phi-\frac{\rho}{b}\sin\phi\Biggr{)}\;\bm{n}\,\Biggr{]},$$

$$\rho =\sqrt{b^2\left(\frac{1+{\tan }^2\varphi }{1+\gamma {\tan }^2\varphi }\right)-\frac{1}{\gamma }{\left(\frac{d}{2}\right)}^2}.$$

$${\mathbf{F}}_0+{\mathbf{F}}_1$$

$$=-\left(\sigma \hat{\mathbf{I}}+\alpha \mathbf{n}\mathbf{n}\right)\cdot {\mathbf{d}}^{\perp },$$

$$F_{\perp }^2+\gamma F_{\parallel }^2=const.,$$

$${\lambda }^3\kappa =-{\lambda }_{\min }^2(\alpha +\sigma )=\mathrm{const}.$$

$$\lambda \kappa =-{\left(\frac{{\lambda }_{\min }}{\lambda }\right)}^2(\alpha +\sigma ).$$

$$Q_{ij}=S\left(n_i{n}_j-\frac{1}{2}{\delta }_{ij}\right),$$

$$\hat{\mathbf{Q}}=\left[\begin{array}{cc}{Q}_{xx}& Q_{xy}\\ Q_{xy}& -Q_{xx}\end{array}\right]=\frac{S}{2}\left[\begin{array}{cc}\cos 2{\theta }_{\mathrm{SF}}& \sin 2{\theta }_{\mathrm{SF}}\\ \sin 2{\theta }_{\mathrm{SF}}& -\cos 2{\theta }_{\mathrm{SF}}\end{array}\right].$$

$${\partial }_t{Q}_{ij}=-\frac{1}{{\xi }_r}\frac{\delta F_{\mathrm{cyto}}}{\delta Q_{ij}},$$

$$\begin{gathered} F_{\mathrm{cyto}}=\frac{K}{2}\int \mathrm{d}A\left[|\nabla \hat{\mathbf{Q}}{|}^2+\frac{1}{{\delta }^2}\mathrm{tr}{\hat{\mathbf{Q}}}^2(\mathrm{tr}{\hat{\mathbf{Q}}}^2-1)\right] \\ +\frac{W}{2}\oint \mathrm{d}s\mathrm{tr}\left[(\hat{\mathbf{Q}}-{\hat{\mathbf{Q}}}_0{)}^2\right]. \end{gathered}$$

$$Q_{0,ij}=T_i{T}_j-\frac{1}{2}{\delta }_{ij},$$

$$K\mathbf{N}\cdot \nabla Q_{ij}-2W\left(Q_{ij}-Q_{0,ij}\right)=0.$$

Full text

I Introduction

Mechanical cues play a vital role in many cellular processes, such as durotaxis Lo2000 ; Sochol2011 , cell-cell communication King2008 , stress-regulated protein expression Sawada2006 or rigidity-dependent stem cell differentiation Engler2006 ; Trappmann2012 . Whereas mechanical forces can directly activate biochemical signaling pathways, via the mechanotransduction machinery Panciera2017 , their effect is often mediated by the cortical cytoskeleton, which, in turn, affects and can be affected by the geometry of the cell envelope.

By adjusting their shape, cells can sense the mechanical properties of their microenvironment and regulate traction forces Bischofs2009 ; Ghibaudo2009 ; Fletcher2010 , with prominent consequences on bio-mechanical processes such as cell division, differentiation, growth, death, nuclear deformation and gene expression Minc2011 ; Versaevel2012 ; Chen1997 ; Jain2013 ; McBeath2004 ; Kilian2010 . On the other hand, the cellular shape itself depends on the mechanical properties of the environment. Experiments on adherent cells have shown that the stiffness of the substrate strongly affects cell morphology Chopra2011 ; Engler2004 and triggers the formation of stress fibers Yeung2005 ; Grinnell2000 . The cell spreading area increases with the substrate stiffness for several cell types, including cardiac myocytes Chopra2011 , myoblasts Engler2004 , endothelial cells and fibroblasts Yeung2005 , and mesenchymal stem cells Zemel2010a .

In our previous work Pomp2018 we have investigated the shape and traction forces of adherent cells Schwarz2013 characterized by a highly anisotropic actin cytoskeleton. Using a contour model of cellular adhesion Bar-Ziv1999 ; Bischofs2008 ; Bischofs2009 ; Schwarz2013 ; Giomi2019 , we demonstrated that the edge of these cells can be accurately approximated by a collection of elliptical arcs obtained from a unique ellipse, whose eccentricity depends on the degree of anisotropy of the contractile stresses arising from the actin cytoskeleton. Furthermore, our model quantitatively predicts how the anisotropy of the actin cytoskeleton determines the magnitudes and directions of traction forces. Both predictions were tested in experiments on highly anisotropic fibroblastoid and epithelioid cells Danen2002 supported by stiff microfabricated elastomeric pillar arrays Tan2003 ; Trichet2012 ; VanHoorn2014 , finding good quantitative agreement.

Whereas these findings shed light on how cytoskeletal anisotropy controls the geometry and forces of adherent cells, they raise questions on how anisotropy emerges from the three-fold interplay between isotropic and directed stresses arising within the cytoskeleton, the shape of the cell envelope and the geometrical constraints introduced by focal adhesions. It is well known that the orientation of the stress fibers in elongated cells strongly correlates with the polarization direction of the cell Vignaud2012 ; Ladoux2016 ; Lam2016 ; Gupta2019 . Consistently, our results indicate that, in highly anisotropic cells, stress fibers align with each other and with the cell’s longitudinal direction (see, e.g., Fig. 1A) Pomp2018 . However, the physical origin of these alignment mechanisms is less clear and inevitably leads to a chicken-and-egg causality dilemma: do the stress fibers align along the cell’s axis or does the cell elongate in the direction of the stress fibers?

The biophysical literature is not scarce of cellular processes that might possibly contribute to alignment of stress fibers with each other and with the cell edge. Mechanical feedback between the cell and the extracellular matrix or substrate, for instance, has been shown to play an important role in the orientation and alignment of stress fibers Deshpande2007 ; Walcott2010 ; Zemel2010a ; Zemel2010b ; Nisenholz2014 . Molecular motors have also been demonstrated to produce an aligning effect on cytoskeletal filaments, both in mesenchymal stem cells Raab2012 and in purified cytoskeletal extracts Schaller2010 , where the observation is further corroborated by agent-based simulations Kraikivski2006 . A similar mechanism has been theoretically proposed for microtubules-kinesin mixtures Swaminathan2009 . Studies in microchambers demonstrated that steric interactions can also drive alignment of actin filaments with each other and with the microchamber walls Soares2011 ; Alvarado2014 ; Deshpande2012 . Theoretical studies have highlighted the importance of the stress fibers’ assembly and dissociation dynamics Deshpande2006 ; Deshpande2007 , the dynamics of focal adhesion complexes Deshpande2008 ; Rens2018 , or both Pathak2008 ; Ronan2014 . Also the geometry of actin nucleation sites has been shown to affect the growth direction of actin filaments, thus promoting alignment Reymann2010 ; Letort2015 . Finally, mechanical coupling between the actin cytoskeleton and the plasma membrane has been theoretically shown to lead to fiber alignment, as bending moments arising in the membrane result into torques that reduce the amount of splay within the filaments Liu2008 . Despite such a wealth of possible mechanisms, it is presently unclear which one or which combination is ultimately responsible for the observed alignment of stress fibers between each other and with the cell’s longitudinal direction. Analogously, it is unclear to what extent these mechanisms are sensitive to the specific mechanical and topographic properties of the environment, although some mechanisms rely on specific environmental conditions that are evidently absent in certain circumstances (e.g., the mechanical feedback between the cell and the substrate discussed in Refs. Thomopoulos2005 ; Deshpande2007 ; Zemel2010b relies on deformations of the substrate and is unlikely to play an important role in experiments performed on stiff micro-pillar arrays).

In this paper we investigate the interplay between the anisotropy of the actin cytoskeleton and the shape of cells adhering to stiff microfabricated elastomeric pillar arrays Tan2003 ; Trichet2012 ; VanHoorn2014 . Rather than pinpointing a specific alignment mechanism, among those reviewed above, we focus on the interplay between cell shape and the spatial organization of the actin cytoskeleton. This is achieved by means of a phenomenological treatment of the stress fiber orientation based on the continuum description of nematic liquid crystals, coupled with a contour model of the cell edge Pomp2018 . The paper is organized as follows: in Sec. II we review in detail our contour model of cells with an anisotropic cytoskeleton Pomp2018 . We study how the anisotropy affects the curvature of and tension in the cell cortex, as well as the forces that the cell exerts on the substrate. In Sec. III we further generalize this approach by studying the mechanical interplay between the orientation of the actin cytoskeleton, modeled as a nematic liquid crystal under confinement of the cell edge, and the shape of the cell, and we compare our results to experimental data on highly anisotropic cells. In both sections we assume that the coordinates of the adhesion sites are constant in time and known. A theoretical description of the dynamics of these adhesion sites, as a result of focal adhesion dynamics, is beyond the scope of this study and can be found, for example, in Refs. Deshpande2008 ; Rens2018 .

II Equilibrium configuration of the cell boundary

Many animal cells spread out after coming into contact with a stiff adhesive surface. They develop transmembrane adhesion receptors at a limited number of adhesion sites that lie mainly along the cell contour (i.e., focal adhesions Burridge1996 ). These cells are then essentially flat and assume a typical shape characterized by arcs which span between the sites of adhesion, while forces are mainly contractile Schwarz2013 . This makes it possible to describe adherent cells as two-dimensional contractile films, whose shape is unambiguously identified by the position $\bm{r}=(x,y)$ of the cell contour Bar-Ziv1999 ; Bischofs2008 ; Bischofs2009 ; Banerjee2013 ; Giomi2013 ; Pomp2018 ; Giomi2019 . Fig. 1B illustrates a schematic representation of the cell (fibroblastoid) in Fig. 1A, where the cell contour consists of a collection of curves, referred to as “cellular arcs”, that connect two consecutive adhesion sites. These arcs are parameterized by the arclength $s$ as curves spanned counterclockwise around the cell, oriented along the tangent unit vector $\bm{T}=\partial_{s}\bm{r}=(\cos\theta,\sin\theta)$, with $\theta=\theta(s)$ the turning angle of the arc with respect to the horizontal axis of the frame of reference. The normal vector $\bm{N}=\partial_{s}\bm{r}^{\perp}=(-\sin\theta,\cos\theta)$, with $\bm{r}^{\perp}=(-y,x)$, is directed toward the interior of the cell. The equation describing the shape of a cellular arc is obtained upon balancing all the conservative and dissipative forces experienced by the cell contour. These are:

[TABLE]

where $t$ is time, $\xi_{t}$ is a (translational) friction coefficient, $\bm{\hat{\Sigma}}_{\rm out}$ and $\bm{\hat{\Sigma}}_{\rm in}$ are the stress tensors on the two sides of the cell boundary and $\bm{F}_{\rm cortex}$ is the stress resultant along the cell contour Schwarz2013 ; Bischofs2008 ; Bischofs2009 ; Banerjee2013 ; Pomp2018 ; Giomi2019 . The temporal evolution of the cell contour is then dictated by a competition between internal and external bulk stresses acting on the cell boundary and the contractile forces arising within the cell cortex. We assume the dynamics of the cell contour to be overdamped.

The stress tensor can be modeled upon taking into account isotropic and directed stresses. The latter are constructed by treating the stress fibers as contractile force dipoles, whose average orientation $\theta_{\rm SF}$ is parallel to the unit vector $\bm{n}=(\cos\theta_{\rm SF},\sin\theta_{\rm SF})$ (see Fig. 1B). This gives rise to an overall contractile bulk stress of the form Pedley1992 ; Simha2002 :

[TABLE]

where $\bm{\hat{I}}$ is the identity matrix, $\sigma>0$ embodies the magnitude of all isotropic stresses (passive and active) experienced by the cell edge and $\alpha>0$ is the magnitude of the directed contractile stresses and is proportional to the local degree of alignment between the stress fibers, in such a way that $\alpha$ is maximal for perfectly aligned fibers, and vanishes if these are randomly oriented. In Sec. III we will explicitly account for the local orientational order of the stress fibers using the language of nematic liquid crystals. The degree of anisotropy of the bulk stress is thus determined by the ratio between the isotropic contractility $\sigma$ and the directed contractility $\alpha$. Finally, the tension within the cell cortex is modeled as $\bm{F}_{\rm cortex}=\lambda\bm{T}$, where line tension $\lambda$ embodies the contractile forces arising from myosin activity in the cell cortex. This quantity varies, in general, along an arc and can be expressed as a function of the arclength $s$. In the presence of anisotropic bulk stresses, in particular, $\lambda(s)$ cannot be constant, as we will see in Sec. II.1. The force balance condition, Eq. (1), becomes then

[TABLE]

II.1 Equilibrium cell shape and line tension

In this section we describe the position of the cell boundary under the assumption that the timescale required for the equilibration of the forces in Eq. (3) is much shorter than the typical timescale of cell migration on the substrate (i.e., minutes). Under this assumption, $\partial_{t}\bm{r}=\bm{0}$ and Eq. (3) can be cast in the form:

[TABLE]

where we have used $\partial_{s}\bm{T}=\kappa\bm{N}$, with $\kappa=\partial_{s}\theta$ the curvature of the cell edge. A special situation is obtained when the cytoskeleton is purely isotropic. In this case, first discussed by Bar-Ziv et al. in the context of cell pearling Bar-Ziv1999 and later expanded by Bischofs and coworkers Bischofs2008 ; Bischofs2009 , $\alpha=0$ and Eq. (4) reduces to

[TABLE]

Solving this equation in the two orthogonal directions $\bm{T}$ and $\bm{N}$ leads to two important insights for isotropic cells. First, in the presence of strictly isotropic bulk forces, the line tension $\lambda$ must be constant along a single arc at equilibrium (i.e., $\partial_{s}\lambda=0$). Second, bulk and peripheral contractility, that pull the cell edge inward and outward along the normal direction, compromise along an arc of constant curvature $\kappa=-\sigma/\lambda$, with the negative sign of $\kappa$ indicating that the arcs are curved inwards. This mechanism is analogous to the Young-Laplace law for fluid interfaces, where the isotropic bulk contractility $\sigma$ plays the role of the pressure difference and the line tension $\lambda$ that of surface tension. The cell edge is then described by a sequence of circular arcs, whose radius $R=1/|\kappa|=\lambda/\sigma$ depends on the local cortical tension $\lambda$ of the arc. The latter depends on the local myosin concentration and, in general, varies from arc to arc. The possibility of a geometry-dependent cortical tension, originating from the elasticity of the transverse fibers, was also explored in Refs. Bischofs2008 ; Bischofs2009 , to account for an observed correlation between the radius and the length of the cellular arcs. Both geometry-dependent and geometry-independent models successfully describe the shape of adherent cells in the presence of strictly isotropic forces. However, as we showed elsewhere Pomp2018 , these isotropic models are not suited for describing cells that develop strong directed forces due to the anisotropic cytoskeleton originating from actin stress fibers Pellegrin2007 ; Burridge2013 .

The presence of an anisotropic cytoskeleton is modeled by a nonzero directed stress, $\alpha>0$. By virtue of Eq. (4), this results into the fact that the cortical tension $\lambda$ is no longer constant along the cellular arcs, since the directed stress has, in general, a non-vanishing tangential component (i.e., $\bm{n}\cdot\bm{T}\neq 0$). In order to highlight the physical mechanisms described by Eq. (4), we introduce a number of simplifications in the remainder of this section. These will be lifted in Sec. III, where we will consider the most general scenario. First, because the orientation of the stress fibers typically varies only slightly along a single arc Pomp2018 , we assume the orientation of the stress fibers, $\theta_{\rm SF}$, to be constant along a single cellular arc, but different from arc to arc. Furthermore, we assume that the contractile bulk stresses $\alpha$ and $\sigma$ are uniform throughout the cell. This allows us to solve Eq. (4) analytically for a given arc. Without loss of generality, we orient the reference frame such that the stress fibers are parallel to the $y-$axis. Thus, $\theta_{\rm SF}=\pi/2$ and $\bm{n}=\bm{\hat{y}}$. Since $\alpha$, $\sigma$ and $\bm{n}$ are constant along an arc, Eq. (4) can be expressed as a total derivative and integrated directly. This yields

[TABLE]

where $\bm{C}_{1}=(C_{1x},C_{1y})$ is an integration constant. Decomposing Eq. (6) into $x-$ and $y-$directions yields

[TABLE]

Next, taking the ratio of Eqs. (7), using $\tan\theta=\mathrm{d}y/\mathrm{d}x$ and integrating, we obtain a general solution for the shape of the cellular arc subject to a non-vanishing isotropic stress (i.e., $\sigma\neq 0$), namely

[TABLE]

where $C_{2}$ is another integration constant and we have set

[TABLE]

Eq. (8) describes an ellipse centered at $(x_{c},y_{c})$ and whose minor and major semi-axis are $a=\sqrt{\gamma C_{2}}$ and $b=\sqrt{C_{2}}$. The constant $\gamma$ quantifies the anisotropy of the bulk stress: $\gamma=1$ corresponds to the isotropic case, whereas $\gamma=0$ represents purely anisotropic bulk contractility. Using again Eqs. (7), we further obtain an expression for the line tension $\lambda$ as a function of $x$ and $y$:

[TABLE]

Using Eqs. (7) and (8), this can be also expressed as a function of the turning angle $\theta$, namely

[TABLE]

This expression highlights the physical meaning of the integration constant $C_{2}$. The right-hand side of Eq. (10) attains its minimal value ($C_{2}$) where $\theta=0$, hence when the tangent vector is perpendicular to the stress fibers (i.e., $\bm{n}\cdot\bm{T}=0$). Thus $C_{2}=\lambda_{\min}^{2}/\sigma^{2}$, where $\lambda_{\rm min}$ is the minimal tension withstood by the cortical actin. Substituting $C_{2}$ in Eq. (10) yields the tension as a function of the turning angle:

[TABLE]

Eq. (11) describes how the line tension deviates from the minimum value $\lambda_{\min}$ for nonzero $\theta$ in order to balance the tangential component of the directed stress [Eq. (4)]. The maximum value of the line tension is found at $\theta=\pi/2$, where the stress fibers are parallel to the arc, and is given by $\lambda_{\max}=\lambda_{\min}/\sqrt{\gamma}$. We note that these variations in line tension along a single arc do not necessarily have to be regulated by the cell. Instead, they could simply be a result of passive mechanical forces in a way very similar to the space-dependent tension in a simple cable hanging under gravity.

Although the minimal line tension $\lambda_{\mathrm{min}}$ could, in principle, be arc-dependent, for example if the cell cortex displays substantial variations in the myosin densities Bischofs2008 , here we approximate $\lambda_{\rm min}$ as a constant. This approximation is motivated by the fact that our previous in vitro observations of anisotropic epithelioid and fibroblastoid cells did not identify a correlation between the arc length and curvature Pomp2018 , which, on the other hand, is expected if $\lambda_{\mathrm{min}}$ was to vary significantly from arc to arc Bischofs2008 . Hence, $\alpha$, $\sigma$ and $\lambda_{\rm\min}$ represent the material parameters of our model.

Substituting $C_{2}$ in Eq. (8) yields an implicit representation of the plane curve approximating individual cellular arcs, namely

[TABLE]

This equation describes an ellipse centered at the point $(x_{c},y_{c})$ and oriented along the $y-$direction, whose minor and major semi-axes are $a=\lambda_{\min}\sqrt{\gamma}/\sigma$ and $b=\lambda_{\min}/\sigma$ respectively (Fig. 2). For arbitrary stress fiber orientation $\theta_{\rm SF}$, Eq. (12) can be straightforwardly generalized in the form:

[TABLE]

Consistently, the major axis of the ellipse is always parallel to the stress fibers, hence tilted by an angle $\theta_{\rm SF}$ with respect to the $x-$axis (Fig. 2). The direct relation between the contractile forces arising from the cytoskeleton and the shape of the cell is highlighted by the dimensionless parameter $\gamma=\sigma/(\sigma+\alpha)$: on the one hand, $\gamma$ defines the anisotropy of the contractile bulk stress, on the other hand it dictates the anisotropy of the cell shape. The coordinates of the center of the ellipse $(x_{c},y_{c})$ and the angular coordinates of the adhesion sites along the ellipse, $\psi_{0}$ and $\psi_{1}$ in Fig. 2, can be calculated using standard algebraic manipulation and are given in Appendix A.

One of the most striking consequences of the elliptical shape of the cellular arcs is that the local curvature is no longer constant, consistent with experimental observations on epithelioid and fibroblastoid cells in Ref. Pomp2018 . This can be calculated from Eq. (12) using $\kappa=(\mathrm{d}^{2}y/\mathrm{d}x^{2})/[1+(\mathrm{d}y/\mathrm{d}x)^{2}]^{\frac{3}{2}}$ and expressed in terms of the turning angle $\theta$ using Eqs. (7) and (11). This yields

[TABLE]

with the negative sign indicating that the arcs are curved inwards. A cellular arc thus attains its maximal (minimal) absolute curvature, where $\theta=0$ ($\theta=\pi/2$) and the stress fibers are parallel (perpendicular) to the arc tangent vector, namely

[TABLE]

Consistent with experimental evidence Pomp2018 , the radius of curvature of arcs perpendicular to stress fibers is smaller than the radius of curvature of arcs parallel to the stress fiber direction. For the experimental methods and a more detailed comparison between theory and experiment, see Ref. Pomp2018 . We stress that, as long as the contractile stresses arising from the actin cytoskeleton are roughly uniform across the cell (i.e., $\alpha$, $\sigma$ and $\lambda_{\mathrm{min}}$ are constant), all cellular arcs are approximated by a unique ellipse (see, e.g., Fig. 3). Thus, regions of the cell edge having higher and lower local curvature correspond to different portions of the same ellipse, depending on the relative orientation of the local tangent vector and the stress fibers.

II.2 Traction forces

With the expressions for shape of the cellular arcs [Eq. (13)] and cortical tension [Eq. (11)] in hand, we can now calculate the traction forces exerted by the cell via the focal adhesions positioned at the end-points of a given cellular arc (Fig. 2). Calling these $\bm{F}_{0}$ and $\bm{F}_{1}$ and recalling the the cell edge is oriented counter-clockwise, we have $\bm{F}_{0}=-\lambda(\theta_{0})\bm{T}(\theta_{0})$ and $\bm{F}_{1}=\lambda(\theta_{1})\bm{T}(\theta_{1})$, where $\theta_{0}$ and $\theta_{1}$ are the turning angles at the end-points of the arc. For practical applications, it is often convenient to express the position of the adhesion sites in terms of their relative distance $\bm{d}=d(\cos\phi,\sin\phi)$ (Fig. 2). This yields

[TABLE]

where $\bm{n}^{\perp}=(\sin\theta_{\mathrm{SF}},-\cos\theta_{\mathrm{SF}})$ and the length scale $\rho$ is defined as

[TABLE]

The total traction force exerted by the cell can be calculated by summing the two forces associated with the arcs joining at a given adhesion site, while taking into account that the the orientation $\bm{n}$ of the stress fibers is generally different from arc to arc.

Another interesting quantity is obtained by adding the forces $\bm{F}_{0}$ and $\bm{F}_{1}$ from the same arc. Although these two forces act on two different adhesion sites, their sum represents the total net force that a single cellular arc exerts on the substrate. This is given by

[TABLE]

where $\bm{d}^{\perp}=d(-\sin\phi,\cos\phi)$ (Fig. 2). Eq. (II.2) is the force resulting from the integrated contractile bulk stress [see Eq. (1)] and is independent of the line tension $\lambda_{\min}$.

II.3 Mechanical invariants

We conclude this section by highlighting two mechanical invariants, local quantities that are constant along a cellular arc, thus showing the intimate relation between the geometry of the cell and the mechanical forces it exerts on the environment. From Eqs. (16) we find

[TABLE]

where $F_{\parallel}$ and $F_{\perp}$ are the components of the force, parallel and perpendicular to $\bm{n}$, at any point along a same cellular arc. Furthermore, by inspection of Eqs. (14) and (11) we observe that

[TABLE]

From this, we find that the normal component of the cortical force, $\lambda\kappa$ [see Eq. (4)], is then given by

[TABLE]

This relation is an analog of the Young-Laplace law for our anisotropic system. In the isotropic limit, $\alpha=0$ and $\lambda_{\min}=\lambda$, thus we recover the standard force-balance expression $\lambda\kappa=-\sigma$. Eq. (21) shows that the normal force $\lambda\kappa$ decreases with increasing line tension $\lambda$, because an increase in line tension is accompanied by an even stronger decrease in the curvature $\kappa$.

III Interplay between orientation of the cytoskeleton and cellular shape

In this Section we generalize our approach by allowing the orientation of the stress fibers to vary along individual cellular arcs. This is achieved by combining the contour model for the cell shape, reviewed in Sec. II, with a continuous phenomenological model of the actin cytoskeleton, rooted into the hydrodynamics of nematic liquid crystals Gennes_book . This setting can account for the mechanical feedback between the orientation of the stress fibers and the cellular shape and allows us to predict both these features starting from the position of the adhesions sites alone. A treatment of the dynamics of focal adhesions is beyond the scope of this paper and can be found elsewhere, e.g., in Refs. Deshpande2008 ; Rens2018 .

As mentioned in the Introduction, experimental observations, by us Pomp2018 and others Vignaud2012 ; Ladoux2016 ; Lam2016 ; Gupta2019 , have indicated that stress fibers tend to align with each other and with the cell’s longitudinal direction. As we discussed, several cellular processes might contribute to these alignment mechanisms, such as mechanical cell-matrix feedback Zemel2010a ; Zemel2010b ; Deshpande2007 ; Walcott2010 ; Nisenholz2014 , motor-mediated alignment Raab2012 ; Schaller2010 ; Kraikivski2006 ; Swaminathan2009 , steric interactions Soares2011 ; Alvarado2014 ; Deshpande2012 , stress fiber formation and dissociation Deshpande2007 ; Deshpande2006 ; Pathak2008 ; Ronan2014 , focal adhesion dynamics Deshpande2008 ; Rens2018 ; Pathak2008 ; Ronan2014 , the geometry of actin nucleation sites Reymann2010 ; Letort2015 , or membrane-mediated alignment Liu2008 , but it is presently unclear which combination of mechanisms is ultimately responsible for the orientational correlation observed in experiments. Our phenomenological description of the actin cytoskeleton allows us to focus on the interplay between cellular shape and the orientation of the stress fibers, without the loss of generality that would inevitably result from selecting a specific alignment mechanism among those listed above.

III.1 Dynamics of the stress fibers

The actin cytoskeleton is modeled as a nematic liquid crystal confined within the cellular contour. This is conveniently represented in terms of the two-dimensional nematic tensor:

[TABLE]

where $\delta_{ij}$ is the Kronecker delta and $S=\sqrt{2\operatorname{tr}\bm{Q}^{2}}$ is the so called nematic order parameter, measuring the amount of local nematic order. In the context of the actin cytoskeleton, a low $S$ value might result from either a low local density of stress fibers, or from a high density of randomly oriented stress fibers. In the standard $\{\bm{\hat{x}},\bm{\hat{y}}\}$ Cartesian basis, Eq. (22) yields

[TABLE]

Just like the dynamics of the cell edge, the dynamics of the cytoskeleton are assumed to be overdamped, thus:

[TABLE]

where $\xi_{r}$ is a rotational friction coefficient, dictating the rate of the relaxational dynamics, and $F_{\rm cyto}$ is the Landau-de Gennes free-energy Gennes_book :

[TABLE]

The first integral in Eq. (25) corresponds to a standard mean-field free-energy, favoring perfect nematic order (i.e., $S=1$), while penalizing gradients in the orientation of the stress fibers and their local alignment. $K$ is the orientational stiffness of the nematic phase in the one elastic constant approximation. The length-scale $\delta$ sets the size of the boundary layer in regions where the order parameter drops to zero to compensate a strong distortion of the nematic director $\bm{n}$, such as in proximity of topological defects. The second integral, which is extended over the cell contour, is the Nobili-Durand anchoring energy Nobili1992 and determines the orientation of the stress fibers along the edge of the cell, with the tensor $\bm{\hat{Q}}_{0}$ representing their preferential orientation. Experimental evidence form our (Fig. 3 and Ref. Pomp2018 ) and other’s work (e.g., Refs. Vignaud2012 ; Ladoux2016 ; Lam2016 ; Gupta2019 ), suggests that, in highly anisotropic cells, peripheral stress fibers are preferentially parallel to the cell edge. The same trend is recovered in experiments with purified actin bundles confined in microchambers Soares2011 ; Alvarado2014 . In the language of Landau-de Gennes theory, this effect can be straightforwardly reproduced by setting

[TABLE]

with $\bm{T}$ the tangent unit vector of the cell edge, under the additional assumption that nematic order is prominent along the cell contour (i.e., $S=1$). The constant $W>0$ measures the strength of this parallel anchoring. To ensure the free-energy to be minimal at equilibrium (i.e., when $\partial_{t}Q_{ij}=0$), we solve Eq. (24) with Neumann boundary conditions:

[TABLE]

III.2 The dynamics of the cell edge

The relaxational dynamics of the cell contour are governed, in our model, by Eq. (3), which is now lifted from the assumption that the orientation $\bm{n}$ of the stress fibers is uniform along individual cellular arcs. Furthermore, the contractile stress given by Eq. (2) can now be generalized as:

[TABLE]

with $\alpha_{0}$ a constant independent on the local order parameter. Comparing Eqs. (2) and (28) yields $\alpha=\alpha_{0}S$, thus the formalism introduced in this Section allows us to explicitly account for the effect of the local orientational order on the amount of contractile stress exerted by the stress fibers.

Next, we decompose Eq. (3) along the normal and tangent directions of the cell contour. Since the cells considered here are pinned at the adhesion sites and the density of actin along the cell contour is assumed to be constant, tangential motion is suppressed, i.e., $\bm{T}\cdot\partial_{t}\bm{r}=0$. Together with Eq. (28) this yields:

[TABLE]

Eq. (29a) describes then the relaxation of tension $\lambda$ within the cell edge, given its shape, whereas Eq. (29b) describes the relaxation of the cellular shape itself. The variations in the cortical tension might result from a regulation of the myosin activity or simply form a passive response of the cortical actin to the tangential stresses.

Integrating Eq. (29a) then yields the cortical tension along an arc:

[TABLE]

where $\bm{\hat{Q}}=\bm{\hat{Q}}(s)$ varies, in general, along an individual cellular arc. The integration constant $\lambda(0)$, which represents the cortical tension at one of the adhesion sites, is related to the minimal tension $\lambda_{\min}$ withstood by the cortical actin which we used, in Sec. II, as material parameter of the problem. In practice, we first calculate $\lambda$ over an entire arc using a arbitrary guess for $\lambda(0)$. Then, we apply a uniform shift in such a way that the minimal $\lambda$ value coincides with $\lambda_{\min}$.

Combining the dynamics of the cell contour and that of the cell bulk, we obtain the following coupled differential equations:

[TABLE]

These are complemented by the condition that $\bm{r}$ is fixed in a number of specific adhesion sites, the boundary condition given by Eq. (27) for the nematic tensor $\bm{\hat{Q}}$ and the requirement that $\min_{s}\lambda(s)=\lambda_{\min}$ on each arc.

IV Numerical results

Eqs. (31) are numerically solved using a finite difference integration scheme with moving boundary. As we detail in Appendix B, the time-integration is performed iteratively using the forward Euler algorithm by alternating updates of the cell contour and of the bulk nematic tensor. This process is iterated until both the cell shape and the orientation reach a steady state.

To highlight the physical meaning of our numerical results, we introduce two dimensionless numbers, namely the contractility number, Co, and the anchoring number, An. Co is defined as the ratio between the typical distance between two adhesion sites $d$ and the major axis of the ellipse approximating the corresponding cellular arc ($b=\lambda_{\mathrm{min}}/\sigma$, see Sec. II.1):

[TABLE]

and provides a dimensionless measure of the cell contraction (thus of the cell average curvature). The anchoring number, on the other hand, is defined as the ratio between a typical length scale $R$ in which the internal cell structure is confined (not necessarily equal to the distance $d$) and the length scale $K/W$, which sets the size of the boundary layer where $\bm{\hat{Q}}$ crosses over from its bulk configuration to $\bm{\hat{Q}}_{0}$:

[TABLE]

This number expresses the ratio between the anchoring energy, which scales as $WR$ [i.e., last term in Eq. (25)], and the bulk energy, which scales as $K$, thus independently on cell size [i.e., first term in Eq. (25)]. Hence, An represents the competition between boundary alignment (with strength $W$) and bulk alignment (strength $K$) within the length scale of the cell $R$. For $\text{\em An}\ll 1$, bulk elasticity dominates over boundary anchoring and the orientation of the stress fibers in the bulk propagates into the boundary, resulting into a uniform orientation throughout the cell and large deviations from parallel anchoring in proximity of the edge. Conversely, for $\text{\em An}\gg 1$, boundary anchoring dominates bulk elasticity and the orientation of the stress fibers along the cell edge propagates into the bulk, leading to non-uniform alignment in the interior of the cell.

To get insight on the effect of the combinations of these dimensionless parameters on the spatial organization of the cell, we first consider the simple case in which the adhesion sites are located at the corners of squares and rectangles (Sec. IV.1). In Sec. IV.2 we consider more realistic adhesion geometries and compare our numerical results with experimental observations on highly anisotropic cells.

IV.1 Rectangular cells

Fig. 4 shows the possible configurations of a model cell whose adhesion sites are located at the vertices of a square, obtained upon varying An and Co, while keeping $\gamma=\sigma/(\sigma+\alpha_{0})$ constant. The thick black line represents the cell boundary, the black lines in the interior of the cells represent the orientation field $\bm{n}$ of the stress fibers and the background color indicates the local nematic order parameter $S$.

As expected, for low Co values (left column), cells with high An exhibit better parallel anchoring than cells with small An values, but lower nematic order parameter in the cell interior. Interestingly, the alignment of stress fibers with the walls in the configuration with large An value (top left) resembles the configurations found by Deshpande et al. Deshpande2006 ; Deshpande2007 , who specifically accounted for the assembly and dissociation dynamics of the stress fibers. More strikingly, the structure reported in the top left of Fig. 4 appears very similar to those found in experiments of actin filaments in cell-sized square microchambers Soares2011 ; Alvarado2014 , simulations of hard rods in quasi-2D confinement Soares2011 , and results based on Frank elasticity Galanis2006 , where the tendency of the nematic director to align along the square edges competes with that of aligning along the diagonal.

For large Co values (right column of Fig. 4), the cell deviates from the square shape. Interestingly, An and Co do not influence the geometry of the cell independently. For constant Co, i.e., for fixed stress fiber contractility, increasing An leads to higher tangential alignment of the stress fibers with the cell edge, thus increasing An decreases the contractile force experienced by the cell edge, which is proportional to $(\bm{n}\cdot\bm{N})^{2}$ [Eq. (31a)].

Finally, we note that all configurations in Fig. 4 display a broken rotational and/or chiral symmetry. For $\text{\em An}=0$ the stress fibers are uniformly oriented, but any direction is equally likely. For non-zero An, the stress fibers tend to align along either of the diagonals (with the same probability) to reduce the amount of distortion. Upon increasing Co, chirality emerges in the cytoskeleton and in the cell contour (see, e.g., the cell in the middle of the right column in Fig. 4). In light of the recent interest in chiral symmetry breaking in single cells Tee2015 and in multicellular environments Duclos2018 , we find it particularly interesting that chiral symmetry breaking can originate from the natural interplay between the orientation of the stress fibers and the shape of the cell.

To conclude this section, we focus on four-sided cells whose adhesion sites are located at the vertices of a rectangle and explore the effect of the cell aspect ratio (i.e., height/width). Fig. 5 displays three configurations with aspect ratio increasing from 1 to 2. Upon increasing the cell aspect ratio, the mean orientation of the stress fibers switches from the diagonal (aspect ratio 1) to longitudinal (aspect ratio 2), along with an increase in the order parameter in the cell bulk, as can be seen in Fig. 5 by the slightly more red-shifted cell interior. This behavior originates from the competition between bulk and boundary effects. Whereas the bulk energy favors longitudinal alignment, as this reduces the amount of bending of the nematic director, the anchoring energy favors alignment along all four edges alike, thus favoring highly bent configurations at the expense of the bulk elastic energy. When the aspect ratio increases, the bending energy of the bulk in the diagonal configuration increases, whereas the longitudinal state only pays the anchoring energy at the short edges, hence independently on the aspect ratio. Therefore, elongating the cell causes the stress fibers to transition from tangential to longitudinal alignment, with a consequent increase of the nematic order parameter. Interestingly, these observations are consistent with the findings of experiments on actin filaments in cell-sized microchambers Soares2011 ; Alvarado2014 .

IV.2 Cells on micropillar arrays

In order to validate our model experimentally, we compare our numerical results with experiments on fibroblastoid and epithelioid cells Danen2002 plated on stiff micropillar arrays Tan2003 ; Trichet2012 ; VanHoorn2014 . The cells are imaged using spinning disk confocal microscopy (see, e.g., Fig. 6A) and the images are then processed in order to detect the orientation of the stress fibers. Upon coarse-graining the local gradients of the image intensity, we reconstruct both the nematic director $\bm{n}$ (black lines, representing the orientation of the stress fibers) and order parameter $S$ (background color, representing the degree of alignment), as visualized in Fig. 6B. See Appendix C for more detail on the experimental methods and image processing.

Consistent with our results on rectangular cells (Fig. 5), the stress fibers align parallel to the cell’s longitudinal direction and perpendicularly to the cell’s shorter edges. Furthermore, the nematic order parameter is close to unity in proximity of the cell contour, indicating strong orientational order near the cell edge, but drops in the interior. This behavior is in part originating from the lower density of stress fibers around the center of mass of the cell, and in part from the presence of $\pm 1/2$ nematic disclinations away from the cell edge. These topological defects naturally arise from the tangential orientation along the boundary. Albeit not uniform throughout the whole cell contour, thus not sufficient to impose hard topological constraints on the configuration of the director in the bulk (i.e., Poincaré-Hopf theorem), this forces a non-zero winding of the stress fibers, which in turn is accommodated via the formation of one or more disinclinations. As a consequence of the concave shape of the cell boundary, these defects have most commonly strength $-1/2$.

To compare our theoretical and experimental results, we extract the locations of the adhesion sites from the experimental data (see Ref. Pomp2018 ) and use them as input parameters for the simulations. In Figs. 6C-E we show results of simulations of the cell in Figs. 6A,B for increasing An values, thus decreasing magnitude of the length scale $K/W$. Here, we take the length scale $R=23.6\>\mu$m to be the square root of the cell area and we use constant values for the ratios $\lambda_{\mathrm{min}}/\sigma=14.7\>\mu$m and $\sigma/(\sigma+\alpha_{0})=0.40$ as found by an analysis of the elliptical shape of this cell Pomp2018 . Fig. 6C shows the results of a simulation where bulk alignment dominates over boundary alignment (An = 0.33, $K/W=71\>\mu$m), resulting in an approximately uniform cytoskeleton oriented along the cell’s longitudinal direction. The nematic order parameter is also approximatively uniform and close to unity. For larger An values (Fig. 6D, An = 1.7 and $K/W=14\>\mu$m), anchoring effects become more prominent, resulting in distortions of the bulk nematic director and the emergence of a $-1/2$ disclination in the bottom left side of the cell. Upon further increasing An (Fig. 6E, An = 8.0 and $K/W=2.9\>\mu$m), the $-1/2$ topological defect moves towards the interior, as a consequence of the increased nematic order along the boundary. This results in a decrease in nematic order parameter in the bulk of the cell, consistent with our experimental data.

A qualitative comparison between our in vitro (Fig. 6B) and in silico cells (Fig. 6E), highlights a number of striking similarities, such as the overall structure of the nematic director, the configuration of the order parameter (large along the cell edge and in the thin neck at the bottom right of the cell) and the occurrence of a $-1/2$ disclination on the bottom-left side. In order to make the comparison quantitative and infer the value of the phenomenological parameters introduced in this Section, we have further analyzed the residual function

[TABLE]

expressing the departure of the predicted configurations of the nematic tensor, $\bm{\hat{Q}}_{\rm sim}$, from the experimental ones, $\bm{\hat{Q}}_{\rm exp}$. The index $i$ identifies a pixel in the cell and $N$ is the total number of pixels common to both numerical and experimental configurations. By construction, $0\leq\Delta^{2}\leq 1$, with 0 representing perfect agreement. Fig. 7 shows a plot of $\Delta^{2}$ versus the anchoring number An for the cell shown in Fig. 6. Consistent with the previous qualitative comparison, the agreement is best at large An values, indicating a substantial preference of the stress fibers for parallel anchoring along the cell edge. For the cell in Fig. 6, $\Delta^{2}$ is minimized for $\text{\em An}=8.0$ ($\Delta^{2}=0.027$), corresponding to a boundary layer $K/W=2.9\>\mu$m. The corresponding numerically calculated configuration is shown in Fig. 6E. However, the flattening of $\Delta^{2}$ for large An values implies that the actual value of An becomes unimportant once the anchoring torques, which determine the tangential alignment of the stress fibers in the cell’s periphery, outcompete the bulk elastic torques. Therefore, we conclude that the cell illustrated in Fig. 6 is best described by $\text{\em An}\gtrsim 5$, corresponding to a boundary layer $K/W\lesssim 5\,\mu$m.

The same analysis presented above has been repeated for five other cells (Fig. 8). The first column shows the raw experimental data, the second column shows the coarse-grained experimental data, and the third column shows the simulations. For these we used the values of $\lambda_{\mathrm{min}}/\sigma$ and $\sigma/(\sigma+\alpha_{0})$ obtained from a previous analysis of the cell morphology Pomp2018 and the An values found by a numerical minimization of $\Delta^{2}$. Also for these cells $\Delta^{2}$ flattens for large An values, and we estimate $\text{\em An}\gtrsim 3$ and $K/W\lesssim 7\>\mu$m. Similar to the cell in Fig. 6, we observe that the overall structure of the nematic director, including the emergence of $-1/2$ topological defects, is captured well by our approach, although there is no precise correspondence between the locations of the defects in the experiments and the theory.

V Discussion and conclusions

In this article we investigated the spatial organization of cells adhering on a rigid substrate at a discrete number of points. Our approach is based on a contour model for the cell shape Bar-Ziv1999 ; Bischofs2008 ; Bischofs2009 ; Schwarz2013 ; Giomi2019 coupled with a continuous phenomenological model for the actin cytoskeleton, inspired by the physics of nematic liquid crystals Gennes_book . This approach can be carried out at various levels of complexity, offering progressively insightful results. Assuming that the orientation of the stress fibers is uniform along individual cellular arcs (but varies from arc to arc), it is possible to achieve a complete analytical description of the geometry of the cell, in which all arcs are approximated by different portions of a unique ellipse. The eccentricity of this ellipse depends on the ratio between the isotropic and directed stresses arising in the actin cytoskeleton, and the orientation of the major axis of this ellipse is parallel to the stress fibers. This method further allows to analytically calculate the traction forces exerted by the cell on the adhesion sites and compare them with traction force microscopy data.

Lifting the assumption that the stress fibers are uniformly oriented along individual cellular arcs allows one to describe the mechanical interplay between cellular shape and the configuration of the actin cytoskeleton. Using numerical simulations and inputs from experiments on fibroblastoid and epithelioid cells plated on stiff micropillar arrays, we identified a feedback mechanism rooted in the competition between the tendency of stress fibers to align uniformly in the bulk of the cell, but tangentially with respect to the cell edge. Our approach enables us to predict both the shape of the cell and the orientation of the stress fibers and can account for the emergence of topological defects and other distinctive morphological features. These predictions are in good agreement with the experimental data and further offer an indirect way to estimate quantities that are generally precluded to direct measurement, such as the cell’s internal stresses and the orientational stiffness of the actin cytoskeleton.

The success of this relatively simple approach is remarkable given the enormous complexity of the cytoskeleton and the many physical, chemical, and biological mechanisms associated with stress fiber dynamics and alignment Zemel2010a ; Zemel2010b ; Deshpande2006 ; Deshpande2007 ; Walcott2010 ; Nisenholz2014 ; Raab2012 ; Schaller2010 ; Kraikivski2006 ; Swaminathan2009 ; Soares2011 ; Alvarado2014 ; Deshpande2012 ; Deshpande2008 ; Rens2018 ; Pathak2008 ; Ronan2014 ; Reymann2010 ; Letort2015 ; Liu2008 . Yet, the agreement between our theoretical and experimental results suggests that, on the scale of the whole cell, the myriad of complex mechanisms that govern the dynamics of the stress fibers in adherent cells can be effectively described in terms of simple entropic mechanisms, as those at the heart of the physics of liquid crystals.

In addition, our analysis demonstrates that chiral symmetry breaking can originate from the natural interplay between the orientation of the stress fibers and the shape of the cell. A more detailed investigation of this mechanism is beyond the scope of this study, but will represent a challenge in the near future with the goal of shedding light on the fascinating examples of chiral symmetry breaking observed both in single cells Tee2015 and tissues Duclos2018 .

In the future, we plan to use our model to investigate the mechanics of cells adhering to micropatterned substrates that impose reproducible cell shapes Thery2010 , with special emphasis to the interplay between cytoskeletal anisotropy and the geometry of the adhesive patches Thery2006a . These systems are not new to theoretical research, but previous studies have focused on either the cytoskeleton Pathak2008 or on cell shape Albert2014 , rather then on their interaction. Furthermore, our framework could be extended to study the role of cytoskeletal anisotropy in cell motility, for instance by taking into account the dynamics of focal adhesions Deshpande2008 ; Rens2018 , biochemical pathways in the actin cytoskeleton Maree2012 , or cellular protrusions and retractions Segerer2015 . Finally, our approach could be extended to computational frameworks such as vertex models, cellular Potts models, or phase field models Albert2016a , and could provide a starting point for exploring the role of anisotropy in multicellular environments such as tissues Eastwood1998 ; VanderSchaft2011 ; Wartlick2011 ; vanOers2014 ; Santos-Oliveira2015 ; Vijayraghavan2016 ; Saw2017 ; Barton2017 .

author contributions

L.G. designed this study with the help of T.S., R.M.H.M. and E.H.J.D. K.S. and L.G. developed the theoretical framework. K.S. and J.E. developed the code, and J.E. performed the simulations. J.E. and W.P. carried out the analysis of the experimental data. K.S. and L.G. wrote the article. All authors commented on the manuscript.

Acknowledgements.

This work was supported by funds from the Netherlands Organisation for Scientific Research (NWO/OCW), as part of the Frontiers of Nanoscience program (L.G.), the Netherlands Organization for Scientific Research (NWO-FOM) within the program on Barriers in the Brain (W.P. and T.S.; No. FOM L1714M), the Netherlands Organization for Scientific Research (NWO-ALW and NWO-ENW) within the Innovational Research Incentives Scheme (R.M.H.M.; Vidi 2010, No. 864.10.009 and Vici 2017, No. 865.17.004), and the Leiden/Huygens fellowship (K.S.).

Appendix A Angular coordinates of the adhesion sites

With reference to the schematic representation of Fig. 2, the coordinates of the center of the ellipse can be expressed as:

[TABLE]

where the length scale $\rho$ is defined in Eq. (17). From Eqs. (35), standard algebraic manipulations allow us to express the angular coordinate $\psi$ of the adhesion points in the frame of the ellipse (Fig. 2a), namely

[TABLE]

Appendix B Numerical methods

The starting point of our numerical simulations are the positions of the adhesion sites. These are directly determined from the experiments by detecting the pillars along the cell contour that are subject to the largest traction forces (see Ref. Pomp2018 for more details). These adhesion sites are fixed during the simulation. Cellular arcs are parameterized in terms of a discrete number of vertices connected by straight edges in a chain-like manner. The bulk of the cell, representing the cytoskeleton, is defined as the region enclosed by the cell edge, and is discretized as a regularly spaced two-dimensional square lattice with $\bm{\hat{Q}}$ defined at every lattice point.

At $t=0$ the cell consists of an irregular polygon enclosing a random configuration of the nematic tensor. Evidently, this bears no resemblance to a real cell, but reduces the risk of a possible bias in the final configuration. The time-integration is performed by alternating updates of the cell contour and of the bulk nematic tensor. This process is iterated until both the cell edge and the cell bulk reach a steady-state. Our code is available upon request.

The configurations of the cytoskeleton in Figs. 4, 5, 6 and 8 have been visualized with Mathematica Version 11.3 (Wolfram Research, Champaign, IL) using the line integral convolution tool.

B.1 Numerical implementation of the cell bulk

Eq. (31b) is numerically solved via a finite-difference scheme. Time integration is performed using the forward Euler method, whereas spatial derivatives are calculated using the centered difference approximation. In order to calculate derivatives at lattice points located in proximity of the edge, we use the boundary conditions, specified in Eq. (27), to express the values of $Q_{xx}$ and $Q_{xy}$ in a number of ghost points located outside the cells. This is conveniently done upon identifying three possible scenarios, illustrated in Fig. 9. 1) There is a single ghost point on the $x-$ or $y-$axis (Fig. 9A). 2) There are two ghost points, one on each axis (Fig. 9B). 3) There are two ghost points on the same axis and possibly a third one on the other axis (Fig. 9C). In the following, we explain how to address each of these cases.

1) Using the centered difference approximation for the first derivative yields the following expression of the nematic tensor at a ghost point located at $(x\pm\Delta x,y)$ or $(x,y\pm\Delta y)$, with $\Delta x=\Delta y$ the lattice spacing:

[TABLE]

The derivative with respect to $x$ in Eq. (37a) can be calculated from Eq. (27), upon taking $\bm{N}=\pm\bm{\hat{x}}$, where the plus (minus) sign correspond to a ghost point located on the left (right) of the central edge point. Thus $\bm{N}\cdot\nabla Q_{ij}=\pm\partial_{x}Q_{ij}$. Analogously, the derivative with respect to $y$ in Eq. (37b), is approximated as $\bm{N}\cdot\nabla Q_{ij}=\pm\partial_{y}Q_{ij}$, where the plus (minus) sign corresponds to a ghost point located below (above) the central edge point. Combining this with Eq. (27), yields:

[TABLE]

The tensor $Q_{0,ij}$ is evaluated via Eq. (26) using the local orientation of the cell edge.

2) If a given lattice point is linked to ghost points in both the $x-$ and $y-$directions, we evaluate equation (38) for both directions independently as explained in the previous paragraph.

3) If a given lattice point is linked to two ghost points in either the $x-$ or $y-$direction, we employ a forward or backward finite difference approximation for the first spatial derivative of $Q_{ij}$ to evaluate $Q_{ij}$ at the ghost points. This yields:

[TABLE]

Finally, if the given lattice point is also linked to a ghost point on the other axis, this is evaluated independently using Eq. (38).

The lattice points enclosed by the cells continuously change during the course of a simulation, as a consequence of the relaxation of the cell shape. In case of inclusion of a new lattice point that was previously located outside the cell, the associated $Q_{xx}$ and $Q_{xy}$ values are generated by averaging over the nearest neighbours (horizontally and vertically, not diagonally) that were inside the cell during the previous time step.

B.2 Numerical implementation of the cell edge

To calculate the line tension $\lambda$, Eq. (30) is discretized as follows:

[TABLE]

where $\lambda_{0}$ is the line tension at the adhesion site at $s=0$ (position $\bm{r}_{0}$) and $\lambda_{k}$ is the line tension at vertex $k$ (position $\bm{r}_{k}$). $N_{\rm arc}$ is the total number of vertices in which cellular arcs are discretized, and $\lambda_{N_{\rm arc}}$ represents the line tension at the other adhesion site. Furthermore, $\Delta s_{n}=|\bm{r}_{n}-\bm{r}_{n-1}|$, $\bm{T}_{n}=(\bm{r}_{n}-\bm{r}_{n-1})/\Delta s_{n}$, $\bm{N}_{n}=\bm{T}_{n}^{\perp}$ and

[TABLE]

with $\bm{\hat{Q}}_{n}$ and $\bm{\hat{Q}}_{n-1}$ the nematic tensor at the vertices $n$ and $n-1$. These are set equal to $\bm{\hat{Q}}$ at the closest bulk lattice point inside the cell among the four points delimiting the unit cell containing the edge vertices $n$ and $n-1$ respectively. If none of these is inside the cell, we set $Q_{xx,n}=Q_{xy,n}=0$. The quantity $\lambda_{0}$ is calculated in such a way that the minimal $\lambda$ value along an arc equates the input parameter $\lambda_{\min}$, representing the minimal tension withstood by the cortical actin.

Next, the position of the vertices $\bm{r}_{k}$, $k=0,1\ldots\,N_{\rm arc}$ is updated upon integrating Eq. (31a) using the forward Euler method. The curvature and normal vector at vertex $k$, $\kappa_{k}$ and $\bm{N}_{k}$, are found by constructing a circle with radius $R$ through vertices $k-1$, $k$, and $k+1$. The vector from vertex $k$ to the center of the circle is then equated to $\pm R\bm{N}_{k}$, with the sign such that $\bm{N}_{k}$ is an inward pointing normal vector, and $\kappa_{k}=\pm 1/R$, with a negative sign for a concave shape and a positive sign for a convex shape. Along each arc, $\bm{r}_{0}$ and $\bm{r}_{N_{\rm arc}}$ represent the positions of the adhesion sites and are kept fixed during simulations.

Appendix C Experimental data analysis

The coarse-grained experimental data displayed in Figs. 6B and 8F-J are obtained using the following four-step procedure.

1) Epithelioid GE11 and fibroblastoid GD25 cells expressing either $\upalpha 5\upbeta 1$ or $\upalpha\rm v\upbeta 3$ (GD$\upbeta$1, GD$\upbeta$3, GE$\upbeta$1 and GE$\upbeta$3) Balcioglu2015 are cultured on microfabricated elastomeric pillar arrays Tan2003 ; Trichet2012 ; VanHoorn2014 . The resulting cells are imaged using spinning disk confocal microscopy and displayed in Figs. 6A and 8A-E. For details, see Ref. Pomp2018 .

2) The locations of the cell interior and the cell edge are found by applying a low-pass filter on the images using Matlab. The interior of the cell is then sampled by overlaying a square lattice of $512\times 512$ pixels ($1\,{\rm pixel}=0.138\times 0.138$ $\mu{\rm m}^{2}$) on the microscope field-of-view (Figs. 6A and 8A-E).

3) For all pixels that are inside the cell, we calculate the nematic tensor using ImageJ with the OrientationJ plugin OrientationJ in the following way. Given the intensity $I(x_{0},y_{0})$ of the image at the point $(x_{0},y_{0})$, we define the symmetric $2\times 2$ matrix $\bm{\hat{J}}=\langle\nabla I\nabla I\rangle$, where $\langle\cdots\rangle=\int w(x,y){\rm d}x\,{\rm d}y\,(\cdots)$ represents a weighted average with $w(x,y)$ a Gaussian with a standard deviation of five pixels ($0.69\>\mu$m) centered at $(x_{0},y_{0})$. The $\bm{\hat{J}}$ matrix can be expressed as:

[TABLE]

where $\Lambda_{\mathrm{max}}$ and $\Lambda_{\mathrm{min}}$ are the largest and smallest eigenvalues of $\bm{\hat{J}}$, $\bm{e}_{\mathrm{min}}$ the eigenvector corresponding to $\Lambda_{\mathrm{min}}$, and $\bm{\hat{I}}$ the two-dimensional identity matrix. The $\bm{\hat{J}}$ matrix is then used to estimate the average stress fiber direction $\bm{u}$:

[TABLE]

Here, the quantity $\bm{\hat{I}}-\langle\bm{u}\bm{u}\rangle$ reflects that the largest gradients in intensity are perpendicular to the orientation of the stress fibers and $\langle|\nabla I|^{2}\rangle=\operatorname{tr}\bm{\hat{J}}=\Lambda_{\max}+\Lambda_{\min}$. Combining Eqs. (42) and (43), we obtain

[TABLE]

Comparing this with the definition of the nematic tensor Gennes_book :

[TABLE]

we find:

[TABLE]

4) The data are further coarse-grained in blocks of $8\times 8$ pixels corresponding to regions of size $1.104\times 1.104\,\mu m^{2}$ in real space. This results in a new $64\times 64$ lattice. The value of the nematic tensor in the new coarse-grained pixels is obtained from an average over those of the original $8\times 8$ pixels located inside the cell. In turn, the coarse-grained pixels are considered inside the cell if more than half of the original pixels were inside the cell.