Response of cells on a dense array of micro-posts

We have analysed the response of cells on a bed of micro-posts idealized as a Winkler foundation using a homeostatic mechanics framework. The framework enables quantitative estimates of the stochastic response of cells along with the coupled analysis of cell spreading, contractility and mechano-sensitivity. In particular the model is shown to accurately predict that: (i) the extent of cell spreading, actin polymerisation as well as the traction forces that cells exert increase with increasing stiffness of the foundation; (ii) the traction forces that cells exert are primarily concentrated along the cell periphery; and (iii) while the total tractions increase with increasing cell area the average tractions are reasonably independent of cell area, i.e. for a given substrate stiffness, the average tractions that are normalized by cell area do not vary strongly with cell size. These results thus suggest that the increased foundation stiffness causes both the cell area and the average tractions that the cells exert to increase through higher levels of stress-fibre polymerization rather than the enhanced total tractions being directly linked through causation to the larger cell areas. A defining feature of the model is that its predictions are statistical in the form of probability distributions of observables such as the traction forces and cell area. In contrast, most existing models present solutions to specific boundary value problems where the cell morphology is imposed a priori. In particular, in line with observations we predict that the diversity of cell shapes, sizes and measured traction forces increase with increasing foundation stiffness. The homeostatic mechanics framework thus suggests that the diversity of observations in in vitro experiments is inherent to the homeostatic equilibrium of cells rather than being a result of experimental errors.

Abstract We have analysed the response of cells on a bed of micro-posts idealized as a Winkler foundation using a homeostatic mechanics framework. The framework enables quantitative estimates of the stochastic response of cells along with the coupled analysis of cell spreading, contractility and mechanosensitivity. In particular the model is shown to accurately predict that: (i) the extent of cell spreading, actin polymerisation as well as the traction forces that cells exert increase with increasing stiffness of the foundation; (ii) the traction forces that cells exert are primarily concentrated along the cell periphery; and (iii) while the total tractions increase with increasing cell area the average tractions are reasonably independent of cell area, i.e. for a given substrate stiffness, the average tractions that are normalized by cell area do not vary strongly with cell size. These results thus suggest that the increased foundation stiffness causes both the cell area and the average tractions that the cells exert to increase through higher levels of stress-fibre polymerization rather than the enhanced total tractions being directly linked through causation to the larger cell areas. A defining feature of the model is that its predictions are statistical in the form of probability distributions of observables such as the traction forces and cell area. In contrast, most existing models present solutions to specific boundary value problems where the cell morphology is imposed a priori. In particular, in line with observations we predict that the diversity of cell shapes, sizes and measured traction forces increase with increasing foundation stiffness. The homeostatic mechanics framework thus suggests that the diversity of observations in in vitro experiments is inherent to the homeostatic equilibrium of cells rather than being a result of experimental errors.

Introduction
Living cells sense and respond to their mechanical, chemical and topological environments. For example, Engler et al. [1] have shown that cell spreading is strongly affected by substrate stiffness with the projected cell area increasing with increasing substrate stiffness. Similar observations showing that cell shape as characterized by aspect ratio is also dependent on substrate stiffness [2] have also been reported. This dependence of direct visual observables to substrate stiffness is linked to the forces that cells exert on the substrate with the total traction forces also increasing with increasing substrate stiffness [3,4]. However, while the statistics of these behaviours are highly reproducible, cells display a fluctuating response that results in a diversity of observables in nominally identical tests. Importantly, this experimental variability is not only a function of the cell type but also a function of the environment with the standard deviation in all the quantities mentioned above decreasing with decreasing substrate stiffness. The variability in direct observables such as cell shape, area and cytoskeletal protein arrangements is also linked to other critical cell functionality. In particular, mechanical, geometric and topological cues direct the differentiation of mesenchymal stem cells (MSCs) [5][6][7]. However, this behaviour is stochastic such that MSCs differentiate mainly but not exclusively into bone cells when cultured on stiff substrates while the probability to differentiate into neuronal cells increases on soft substrates [5]. Thus, the observed responses of cells are always characterised in terms of statistics rather than unique outcomes. A mechanistic understanding of this stochastic behaviour of cells will have farreaching implications in aiding the interpretation of a wide range of cell functionalities and responses.
Cytoskeletal tension plays an essential role in the way cells sense and respond to their environment with the above discussed sensitivities suppressed when inhibitors such as cytochalasin D are added to diminish myosin contractility [5,8,9]. This understanding has spurred extensive experimental activity to measure the traction forces that cells exert. The main approaches include traction force microscopy [10,11] and using deformable micro-post arrays [12][13][14][15]. The key conclusions from these studies are: (i) the so-called total traction forces rise with increasing substrate stiffness; and (ii) cells spread to a greater extent on stiffer substrates. In fact, both these processes are coupled in the sense that cells undergo remodelling and reorganize their cytoskeleton in response to their mechanical environment, with the cytoskeletal structure controlling both cell shape/area as well as the traction forces that are generated.
A number of modelling approaches have been developed to model the mechano-sensitive response of cells with the aim of using models to try and interpret the traction force measurements. Early attempts included modelling the cytoskeleton as an interlinked structure of passive filaments [16] or modelling contractility as a thermal contraction [17]. Deshpande et al. [18] introduced a bio-chemo-mechanical model for cell contractility that subsequently has been adapted in a number of studies [19][20][21]. These models rationalised a number of the observations of cells on elastic substrates but inherently solved a boundary value problem such that the size and shape of the cell was specified a priori. In reality, contractility and the shape that cells adopt by spreading are interlinked with both, in turn, depending on environmental parameters such as substrate stiffness and ligand density [1]. Comprehensive intricate coupling among cell configuration, traction forces and the elastic environment is missing in all these modelling approaches. Moreover, they necessarily give only deterministic predictions for the response of cells while the measured responses are clearly stochastic.
The key assumption in all the models mentioned above is that they presume, either implicitly or explicitly, that a system comprising the cell and the substrate attains equilibrium at its minimum value of Gibbs free-energy. However, such a system that excludes the nutrient bath surrounding the cell is an open system with the cell exchanging (chemical) species with the nutrient bath (Fig. 1a). In fact, thermodynamic equilibrium of this open system is achieved when the chemical potentials of all mobile species within the cell and nutrient bath equalise but living cells never achieve such an equilibrium state (e.g. all living cells maintain a resting potential between the cell and the surrounding nutrient bath by actively regulating the concentration of various ions within the cell). Hence, cells are inherently in a non-equilibrium state from a conventional thermodynamic perspective. Shishvan et al. [22] have recently proposed a homeostatic ensemble to describe the dynamic equilibrium of cells in an attempt to resolve this critical limitation in all previous modelling approaches. This framework therefore inherently allows us to make quantitative predictions of the stochastic response of cells with cell shape and size no longer imposed via a boundary value problem but rather the distribution of configurations that cells adopt being an outcome of the homeostatic equilibrium in the given environment.
The homeostatic statistical mechanics description for cells [22] has already been shown to successfully capture a range of observations for smooth muscle cells seeded on elastic substrates [22,23] and for myofibroblasts seeded on substrates micropatterned with stripes of fibronectin [24,25] as well as for the differentiation of hMSCs in response to a range of environmental cues including stiffness of substrates and sizes of adhesive islands [26]. These give us confidence in utilizing the homeostatic mechanics framework to investigate the response of cells on a dense array of micro-posts.

Cells on a dense array of micro-posts
Here we consider the problem of a cell on a bed of micro-posts arranged in a cubic pattern with spacing L, as sketched in Fig. 1a. The posts have a height H and Fig. 1 a Sketch of a cell on a bed of micro-posts. b Sketch of the cell on a Winkler foundation used to approximate the bed of micro-posts. In a and b the nutrient bath has been illustrated and a small selection of the species being exchanged between the cell and the bath are labelled. c The two-dimensional (2D) approximation of a cell on the Winkler foundation analysed in this study in line with most experimental studies [12][13][14][15] are assumed to have a circular cross-section of radius a, as sketched in Fig. 1a. The bending stiffness S of the posts is then where E is the Young's modulus of micro-posts. If a representative spread dimension R of the cell in the x 1 À x 2 plane is much greater than the post spacing, i.e. R ) L, then it suffices to approximate the bed of posts as an effective Winkler foundation with smeared-out isotropic stiffness j ¼ S=L 2 in the x 1 À x 2 plane (Fig. 1b). In this study we shall analyse the behaviour of smooth muscle cells (SMCs) on a dense array of micro-posts using this Winkler approximation. We emphasize that the Winkler-foundation approximation implies that the precise post geometry does not enter the analysis explicitly. Rather the analysis only depends on the Winkler stiffness j. The outline of the paper is as follows. We first briefly review the model of Vigliotti et al. [27] as modified in [22] to calculate the free-energy of a system comprising the cell and foundation. Living cells, however, do not attain thermodynamic equilibrium and thus do not equilibrate at a minimum freeenergy state. Next, we briefly describe the homeostatic mechanics framework [22] to account for the dynamical equilibrium that living cells attain. Finally, we present predictions using the homeostatic mechanics framework for the responses of cells on arrays of micro-posts in an attempt to identify correlations or the lack of them in the dependences of cell area, cell shape and traction forces on the stiffness of the micropost arrays.
2 Free-energy of a cell on a dense array of posts Consider a cell in a nutrient bath on a bed of microposts (approximated as a Winkler foundation), as sketched in Fig. 1b. We define a system comprising the cell and foundation but absent the nutrient bath which is considered as the environment. This of course has the advantage that it is not then required to model the nutrient bath but rather only model the parts of the experimental setup on which observations are being directly conducted. The state of the system changes as the cell moves, spreads and changes shape on the foundation and here we shall give a prescription to calculate the free-energy of the cell in a specific configuration ðjÞ where the connections of material points on the cell membrane to the foundation surface are specified (with the remainder of the cell surface and foundation assumed to be traction-free). In broad terms the adhesion of the cell membrane to the foundation specifies the shape of the cell and subsequently we shall refer to each such configuration as a morphological microstate ðjÞ of the system.

The model for the Gibbs free-energy
With the system comprising of the cell and the foundation within a constant temperature and pressure nutrient bath, the Gibbs free-energy G ðjÞ of the system in morphological microstate ðjÞ is given by where f is the specific Helmholtz free-energy of the cell and w the strain energy density of the foundation. We emphasize that the analysis presented here is for the system under atmospheric pressure conditions and thus without loss of generality we set p ¼ 0 (i.e. use gauge pressure). Thus, a pressure term does not appear in (2). The equilibrium free-energy G ðjÞ is then the value of G ðjÞ at dG ðjÞ ¼ 0. Here, we briefly describe the model for the calculation of G ðjÞ . In the following, for the sake of notational brevity, we shall drop the superscript ðjÞ that denotes the morphological microstate as the entire discussion refers to a single morphological microstate.
With the cell exerting a spatial distribution of tractions T 1 ðx i Þ and T 2 ðx i Þ on the foundation, the strain energy density w is given by The total foundation strain energy where the spatial integration is carried out over the current cell area A. We now proceed to summarise the model of Vigliotti et al. [27] for the cell as modified in [22] for a non-dilute concentration of stress-fibres.
Here, we restrict ourselves to a two-dimensional (2D) approximation of the cell (Fig. 1c) such that the cell lies in the x 1 À x 2 plane with no variation of properties in the through thickness x 3 -direction. The Vigliotti et al. [27] model assumes only two elements within the cell: (i) a passive elastic contribution from elements such as the cell membrane, intermediate filaments and microtubules and (ii) contractile acto-myosin stress-fibres that are modelled explicitly. Consider a 2D cell of thickness b 0 and volume V 0 in its elastic resting state. The representative volume element (RVE) of the stress-fibres within the cell in this resting configuration is assumed to be a cylinder of volume V R ¼ pb 0 where ' 0 is the length of a stress-fibre functional unit in its groundstate and n R is the number of these ground-state functional units within this reference RVE. The total number of functional unit packets within the cell is N T 0 and we introduce N 0 ¼ N T 0 V R =V 0 as the average number of functional unit packets available per RVE; N 0 shall serve as a useful normalisation parameter. The state of the stress-fibres at location x i within the cell is described by their angular concentration gðx i ; /Þ, and there are nðx i ; /Þ functional units in series along the length of each stress-fiber in the RVE. Here, / is the angle with respect to the x 1 direction (Fig. 1c). Vigliotti et al. [27] argue that an applied stretch is shared equally among all subunits, so that the strain within each functional unit,ẽ nom , is initially equal to the nominal strain e nom ðx i ; /Þ in direction /. Subsequent addition or removal of subunits modifies the subunit stretch proportionally so that, at steadystate, the number n ss of functional units within the stress fibers is given bŷ whereẽ ss nom is the strain at steady-state within a functional unit of the stress-fibres. It now remains to specify the steady-state angular concentration of the stress-fibres. The chemical potential of the functional units within the stress-fibres is given by [22] v where k and T are Boltzmann constant and temperature, respectively, and the normalized concentration of the unbound stress fiber proteins is given bŷ N u N u =N 0 . Withĝ gn R =N 0 , the maximum normalised value ofĝ, denoted byĝ max , corresponds to full occupancy of all available sites for stress-fibres. Here,N L is the number of lattice sites available to unbound proteins. The enthalpy l ss b of n R bound functional units at steady-state is given in terms of the isometric stress-fibre stress r max and the internal energy l b0 as where X is the volume of n R functional units. By contrast, the chemical potential of the unbound proteins is independent of stress and given in terms of the internal energy l u as For a fixed configuration of the 2D cell (i.e. a fixed strain distribution e nom ðx i ; /Þ), the contribution to the specific Helmholtz free-energy of the cell, f , from the stress-fibre cytoskeleton follows as where q 0 N 0 =V R is the number of protein packets per unit reference volume available to form functional units in the cell. However, we cannot yet evaluate f cyto asN u ðx i Þ andĝðx i ; /Þ are unknown. These will follow from the chemical equilibrium of the cell as will be discussed in Sect. 2.2. The total stress R ij within the cell includes contributions from the passive elasticity provided mainly by the intermediate filaments of the cytoskeleton attached to the nuclear and plasma membranes and the microtubules as well as the active contractile stresses of the stress-fibres. Thus, in line with the existing literature on active stress mechanics, we write the total Cauchy stress in an additive decomposition as where r ij and r p ij are the active and passive Cauchy stresses, respectively. In the 2D setting with the cell lying in the x 1 À x 2 plane, the active stress is given in terms of the volume fraction F 0 of the stress-fibre proteins as where / Ã is the angle of the stress-fibre measured with respect to x i and is related to / by the rotation with respect to the undeformed configuration. The passive elasticity in the 2D setting is given by a 2D specialization of the Ogden [28] hyperelastic model as derived in [22]. The strain energy density function of this 2D Ogden model is where k I and k II are the principal stretches, l and K are the shear modulus and in-plane bulk modulus, respectively, and m is a material constant governing the nonlinearity of the deviatoric elastic response. Moreover, since the cell is assumed to be incompressible we set the principal stretch in the x 3 -direction k III ¼ 1=ðk I k II Þ. The (passive) Cauchy stress then follows as r p ij p ðkÞ j ¼ r p k p ðkÞ i in terms of the principal (passive) Cauchy stresses r p k k k oU elas =ok k and the unit vectors p ðkÞ j ðk ¼ I; IIÞ in the principal directions. The total specific Helmholtz free-energy of the cell is We emphasize here that the constitutive formulation differs considerably from the formulations employed for the inelastic deformation of non-active materials. In such non-active materials, a multiplicative decomposition of the deformation gradient in terms of the elastic and inelastic parts is assumed. By contrast, for the active behaviour of cells we assume that there is no inelastic deformation with the total deformation gradient equal to the elastic deformation gradient. However, an active stress resulting from the stress fibres is assumed to act in parallel to the passive elastic stress.

Equilibrium of the morphological microstate
Shishvan et al. [22] have shown that equilibrium of a morphological microstate reduces to two conditions: the chemical potentials of bound and unbound stress-fibre proteins are equal throughout the cell. The condition v u ¼ v b implies thatĝðx i ; /Þ is given in terms ofN u bŷ andN u follows from the conservation of stress-fibre proteins throughout the cell, viz.
KnowingN u andĝðx i ; /Þ, the stress R ij can now be evaluated and these stresses within the system (i.e. cell and foundation) need to satisfy mechanical equilibrium, i.e. R ij;j ¼ 0. In this case, the mechanical equilibrium condition is readily satisfied as the stress field R ij within the cell is equilibrated by a traction field T i exerted by the foundation on the cell such that bR ij;j ¼ ÀT i , where bðx i Þ is the thickness of the cell in the current configuration. Tractions calculated in this manner are then directly used in Eqs. (3) and (4) to obtain the foundation elastic energy.
The equilibrium value of G denoted by G is then and F foun is specified by Eq. (4). Here, v u is given by Eq. (8) with the equilibrium value ofN u obtained from Eq. (14). For the purposes of further discussion, we define the equilibrium value F cyto q 0 V 0 v u as the cytoskeletal free-energy of the cell and F passive R V cell U elas dV as the passive elastic energy of the cell. The free-energy G can be decomposed as is independent of the morphological microstate. It is thus natural to subtract out ! 0 and define a normalised free-energy aŝ where G S is the equilibrium free-energy of a freestanding cell (i.e. a cell in suspension with tractionfree surfaces). Analogously, we define the normalised passive and cytoskeletal free-energies of the cell aŝ and respectively. We shall present all energies in these normalised forms.

Model parameters
The simulations were performed with the parameters taken from Shishvan et al. [22] and are representative of smooth muscle cells (SMCs). All simulations are reported at a reference thermodynamic temperature The passive elastic parameters of the cell are taken to be l ¼ 1:67 kPa, K ¼ 35 kPa and m ¼ 6. For SMCs, the maximum contractile stress r max ¼ 240 kPa consistent with a wide range of measurements on muscle fibres [29] and the density of stress-fibre proteins was taken as q 0 ¼ 3 Â 10 6 lm À3 with the volume fraction of stress-fibre proteins F 0 ¼ 0:032. Following Vigliotti et al. [27], we assume that the steady-state functional unit straiñ e ss nom ¼ 0:35 with l b0 À l u ¼ 2:3kT 0 and X ¼ 10 À7:1 lm 3 . The maximum angular stress fibre concentrations was set to beĝ max ¼ 1 based on the assumption that the local density of bound stress-fibre proteins cannot exceed q 0 . All results are presented for a cell that is assumed to be circular with a radius R 0 and thickness b 0 in its undeformed state with b 0 =R 0 ¼ 0:2. Since results are presented with energies normalised as per (16)- (18), the quantities l u ; n R and N L do not need to be explicitly specified. With these parameters, the equilibrium free-energy of the freestanding cell is ðG S À ! 0 Þ=ðV 0 kT 0 Þ % À5:6 Â 10 6 lm À3 where V 0 ¼ pR 2 0 b 0 is the cell volume. In this free-standing state the equilibrium configuration of the cell is a spatially uniform cylinder with radius $ 0:92R 0 .

Competition between elastic and cytoskeletal free-energy
One of the key features of the constitutive model of Vigliotti et al. [27] is the competition between the elastic free-energy F passive and the cytoskeletal freeenergy F cyto of the cell that sets the free-energy of the cell. To illustrate this competition, we consider the highly simplified problem of a circular cell on a rigid foundation and constrain ourselves to morphological microstates wherein the strain distribution within the cell is spatially uniform. We emphasize here that this is an unrealistic restriction of the phase-space of morphological microstates that the cell will attain and this restriction is relaxed when the homeostatic mechanics is introduced in Sect. 3. However, for the purposes of illustrating the basic physics of the freeenergy model we present this restrictive analysis here in which a morphological microstate is described by one scalar variable, e.g. the area A of the cell. For the case of a cell on a rigid foundation, there is no contribution to the Gibbs free-energy of the system from the foundation and so G ¼ F cell . The normalised free-energy of the systemĜ is plotted in Fig. 2a as a function of the normalised cell areaÂ A=A 0 , where A 0 is the area of the undeformed cell. There is a clear minimum ofĜ atÂ opt % 1:44. To understand this minimum, the variations of the free-energiesF passive andF cyto withÂ are also included in Fig. 2a. The elastic energy increases with increasing jÂ À 1j as the cell is strained away from its undeformed configuration. By contrast,F cyto decreases monotonically with increasingÂ. This competition betweenF passive and F cyto with varyingÂ gives rise to the minimum inĜ. Assuming that the cell is evolving towards a minimum free-energy configuration (this is clearly not true and we shall show in Sect. 3 that low Gibbs free-energy states are just more probable), we can say that the stress-fibre cytoskeleton drives cell spreading. This is consistent with a large number of observations [8,9] that indicate that inhibiting stress-fibres via reagents such as cytochalasin D and blebbistatin reduces cell spreading. In fact, the model predicts reduced spreading as the available stress-fibre proteins are reduced and this is illustrated in Fig. 2b where we show the dependence ofÂ opt on q 0 keeping all other parameters fixed. With decreasing q 0 the areaÂ opt at which the Gibbs free-energy is minimized becomes closer to the undeformed state withÂ opt ¼ 1.
At face value, stress-fibres driving cell spreading is rather counter-intuitive as stress-fibres exert contractile forces and hence one would expect them to contract the cell rather than promote spreading. To understand this apparent contradiction, recall that the number of functional units in the bound state increases with increasing strain as quantified in Eq. (5). This decreases the number of unbound stress-fibre proteinŝ N u that in turn decreases v u and therefore reduceŝ F cyto . Another way to view this is to recall that the enthalpy of functional units in the bound state is lower than their corresponding enthalpy in the unbound state due to the tensile stress r max within the stress-fibres; see Eq. (7). Chemical equilibrium dictates that all stress-fibres proteins are at equal chemical potentials. This immediately implies that the cytoskeletal freeenergy decreases with increasing strain as the bound protein numbers rise with increasing strain. Thus, it is the formation of stress-fibres with tensile stresses that tends to reduce the Gibbs free-energy of the cell and drives cell spreading. Adding reagents such as cytochalasin D that inhibit the formation of stress-fibres will have the effect of diminishing the reduction inF cyto with increasing strain and therefore tend to reduce cell spreading. We emphasize that cell spreading requires kinetic processes such as polymerization of meshwork actin along the cell periphery and in lamellipodia. These processes are not accounted for here but rather we argue that the overall driving force for spreading is the reduction in the Gibbs free-energy of the system and the availability of kinetic pathways is the means of achieving this reduction in the Gibbs free-energy.

Homeostatic mechanics for cells
The system modelled above comprised the cell and the foundation but excluded the nutrient bath that not only maintains the system at a constant temperature and pressure but also furnishes the cell with nutrients. While cells are alive they maintain a resting potential between themselves and the surrounding nutrient bath by actively regulating the concentration of various ions within the cell [30] via a very large number of complex inter-linked metabolic reactions such as (but not restricted to) ion-pumps, osmosis, diffusion and cytoskeletal reactions. Hence, from a conventional thermodynamics perspective cells are inherently in a non-equilibrium state [31]. Remarkably, the metabolic processes maintain the concentrations of all species within the cell to be very nearly constant over a Fig. 2 a The normalised Gibbs free-energyĜ, passive elastic energyF passive and the cytoskeletal free-energyF cyto as a function of the normalised areaÂ of the cell (using the reference properties). The minimum value ofĜ corresponding to an areâ A opt has been labelled. b The dependence ofÂ opt to the number of protein packets q 0 per unit volume available to form functional units in the cell. In both a and b, the morphological microstates are constrained to comprise only uniformly strained circular cells sustained period of time (e.g. the interphase period of the cell cycle) and this phenomenon is known as cellular homeostasis. Shishvan et al. [22] developed a new thermodynamic framework whereby they formally defined the (dynamic) homeostatic equilibrium. We shall use this framework to analyse the distribution of states that cells assume on elastic Winkler foundations that idealise the array of micro-posts. In this section, we briefly review this homeostatic mechanics framework for cell mechanics with readers referred to [22] for details.

Brief overview of the homeostatic ensemble
The homeostatic statistical mechanics description for cells is applicable over a timescale from a few hours to a few days covering the interphase period of the cell cycle when the cell remains as a single undivided entity. Controlling only macro variables (i.e. macrostate) such as the temperature, pressure and nutrient concentrations in the nutrient bath results in inherent uncertainty (referred to here as missing information) in micro variables (i.e. microstates) of the system. This includes a level of unpredictability in homeostatic process variables, such as the spatio-temporal distribution of chemical species, that is linked to Brownian motion and the complex feedback loops in the homeostatic processes. Thus, this system not only includes the usual lack of precise information on the positions and velocities of individual molecules associated with the thermodynamic temperature, but also an uncertainty in cell shape resulting from the homeostatic processes not being precisely regulated. The consequent entropy production forms the basis of this new statistical mechanics framework motivated by the following two levels of microstates: (i) Molecular microstates Each molecular microstate has a specific configuration (position and momentum) of all the molecules within the system. (ii) Morphological microstates (Fig. 1b) Each morphological microstate is specified by the mapping (connection) of material points on the cell membrane to material points on the foundation. In broad terms, a morphological microstate specifies the shape and size of the cell.
In the homeostatic state, the system is in (dynamic) equilibrium with no net change in the internal state of the system but with a net flux of species between the system and nutrient bath (e.g. there is an overall flux of glucose into the cell while the net flow of carbon dioxide is in the opposite direction). Shishvan et al. [22] identified this (dynamic) equilibrium state by entropy maximisation. Thus, subsequently, we shall simply refer to this state as an equilibrium state to emphasise that it is a stationary macrostate of the system inferred via entropy maximisation as in a conventional equilibrium analysis. The total entropy of the system is written in terms of the conditional probability P ðijjÞ of the molecular microstate ðiÞ given the morphological microstate ðjÞ and the probability P ðjÞ of morphological microstate ðjÞ as In Eq. (19), I ðjÞ M À P i2j P ðijjÞ ln P ðijjÞ and I C À P j P ðjÞ ln P ðjÞ are the entropies of molecular microstates in morphological microstate ðjÞ and the morphological microstates, respectively. Equilibrium then corresponds to molecular and morphological macrostates that maximise I T subject to the appropriate constraints. By identifying the fact that the evolution of the molecular and morphological macrostates is temporally decoupled, Shishvan et al. [22] showed that Eq. (19) can be maximised by independently maximising I ðjÞ M at the smaller timescales to determine the equilibrium distribution of molecular microstates (i.e. molecular macrostate) for a given morphological microstate and then maximising I C at the larger timescale to determine the equilibrium distribution of the morphological microstates (i.e. morphological macrostate).
Over the (short) timescale on the order of seconds, the only known constraint on the system is that it is maintained at a constant temperature, pressure and strain distribution. The equilibrium of a given morphological microstate ðjÞ obtained by maximising I ðjÞ M corresponds to molecular arrangements that minimise the Gibbs free-energy with G ðjÞ and S ðjÞ M being the equilibrium values of the Gibbs free-energy and entropy of morphological microstate ðjÞ, respectively. Since the connection between the cell and the foundation is fixed for a given morphological microstate, determination of G ðjÞ is a standard boundary value problem as described in Sect. 2. Over the (long) timescale on the order of hours, the equilibrium distribution P ðjÞ eq is determined by maximising I C but now with the additional constraint that the cell is maintained in its homeostatic state. For the case of a cell on an elastic foundation in a constant temperature and pressure nutrient bath, the homeostatic constraint translates to the fact that the average Gibbs freeenergy of the system, over all the morphological microstates it assumes, is equal to the equilibrium Gibbs free-energy G S of an isolated cell in suspension (free-standing cell), i.e. the homeostatic processes maintain the average biochemical state of the system equal to that of the cell in suspension. In deriving this result, Shishvan et al. [22] did not consider every individual homeostatic process but rather just used the coarse-grained outcome of the homeostatic processes, viz. over the homeostatic state, the average numbers of all species within the cell are maintained at a fixed value independent of the extracellular environment. The application of this coarse-grained constraint is the key element of the homeostatic mechanics framework with the morphological entropy I C parameterising the information lost by not modelling all the variables associated with the homeostatic processes.
The maximisation of I C while enforcing P j P ðjÞ G ðjÞ ¼ G S gives the homeostatic equilibrium state such that where Z P j exp ÀfG ðjÞ À Á is the partition function of the morphological microstates and the distribution parameter f follows from homeostatic constraint The collection of all possible morphological microstates that the system assumes while maintaining its homeostatic equilibrium state is referred to as the homeostatic ensemble. The homeostatic ensemble can therefore be viewed as a large collection of copies of the system, each in one of the equilibrium morphological microstates. While G ðjÞ of each copy is not known exactly, the copies are distributed in the ensemble as per the exponential distribution P ðjÞ eq with the distribution parameter f. A crucial difference between the canonical and homeostatic ensembles is that unlike T in the canonical ensemble, f is not a property of the nutrient bath but rather set by the homeostatic state that the system attains. Of course, from Eq. (20) we see that morphological microstates with lower free-energy are more probable (i.e. some sort of justification for the minimum free-energy discussion presented in Sect. 2.4) but of course the homeostatic ensemble will also contain morphological microstates with much higher free-energies, depending on the distribution parameter f.
The equilibrium morphological entropy S C ¼ À P j P ðjÞ eq ln P ðjÞ eq (i.e. maximum value of I C ) is related to f via the conjugate relation oS C =oG S ¼ f. Thus, analogous to 1=T that quantifies the increase in the uncertainty of the molecular microstates (i.e. molecular entropy S ðjÞ M ) with average enthalpy, f specifies the increase in the uncertainty of the morphological microstates (i.e. morphological entropy S C ) with the average Gibbs free-energy. We therefore refer to 1=f as the homeostatic temperature with the understanding that it quantifies the fluctuations on a timescale much slower than that characterised by T.

Numerical procedure to sample the homeostatic ensemble
We employ Markov Chain Monte Carlo (MCMC) to construct a Markov chain that is representative of the homeostatic ensemble. This involves three steps: (i) a discretisation scheme to represent a morphological microstate ðjÞ, (ii) calculation of G ðjÞ for a given morphological microstate ðjÞ and (iii) constructing the Markov chain comprising these morphological microstates. Here, we briefly describe the procedure which was programmed in MATLAB with readers referred to [22] for further details.
In the general setting of a three-dimensional (3D) cell, a morphological microstate is defined by the connection of material points on the cell membrane to the surface of the foundation. In the 2D context, this reduces to specifying the connection of all material points of the cell to the foundation, i.e. a displacement field u  M). In all the numerical results presented here we employ M ¼ 32 with 4 Â 4 weights governing the displacements in the x 1 and x 2 directions, respectively. The NURBS employ third order base functions for both the x 1 and x 2 directions and the knots vector included two nodes each with multiplicity three, located at the extrema of the interval. We emphasise here that this choice of representing the morphological microstates imposes restrictions on the morphological microstates that will be considered. Therefore, the choice of the discretisation used to represent u ðjÞ i needs to be chosen so as to be able to represent the microstates we wish to sample, e.g. the choice can be based on the minimum width of a filopodium one expects for the given cell type. Given u ðjÞ i we can calculate G ðjÞ using the model described in Sect. 2.
We construct, via MCMC, a Markov chain that serves as a sample of the homeostatic ensemble. This is done using the Metropolis [32] algorithm in an iterative manner using the following procedure: (i) Assume a value of f and use the undeformed cell configuration as the starting configuration and label it as morphological microstate j ¼ 0 with equilibrium free-energy G ð0Þ calculated as described above.
(ii) Randomly pick two of the M weights U ðjÞ L and perturb them by two independent random numbers picked from a uniform distribution over the interval ½ÀD D.
(iii) Compute the new free-energy G ðjÞ of this perturbed state and thereby the change in free-energy DG ¼ G ðjÞ À G ðjÀ1Þ . (iv) Use the Metropolis criterion to accept this perturbed state or not, i.e.
(a) if DG 0, accept the perturbed state; (b) if DG [ 0, compute P acc ¼ expðÀfDGÞ and accept the perturbed state if P acc [ R, where R is a random number drawn from a uniform distribution over ½0 1: (v) If the perturbed state is accepted add it to the list of samples as a new morphological microstate else restore the configuration prior to step (ii) in the sample list and return to step (ii). (vi) Keep repeating this procedure until a converged distribution is obtained. Here, we typically use the criterion that the average of G ðjÞ within the generated sample list (labelled hG ðjÞ i) changes by less than 1% over 100,000 steps of the Markov chain. Typical Markov chains comprised in excess of 4 million samples. (vii) If hG ðjÞ i is within AE2% of G S we will accept this distribution else we will modify f and repeat from step (i).

Predictions of the response of cells on a dense array of micro-posts
We present results for the response of cells on a bed of micro-posts approximated as a Winkler foundation with a normalised stiffnessĵ jR 0 =l (see Sect. 2.3 for details of model parameters). The predictions are representative of experiments [12][13][14][15] where these is a low seeding density of cells and that there is no cellcell interaction. Therefore, the simulations are performed for single cells. Moreover, for the range of Winkler stiffnessesĵ considered here we have confirmed that the cell tractions remain sufficiently small so that the two critical assumptions used within the Winkler foundation stiffness derivation (1), viz.
(i) there is no post-post interaction and (ii) the post deflections are sufficiently small that linear beam theory suffices, remain valid. The aim here is to demonstrate the application of the model in the general sense rather than aim to model a specific cell type or experiments. Hence, in the following we restrict ourselves to qualitative comparisons with existing experimental data. Predictions of the probability density functions pðĜÞ / wðĜÞ expðÀfĜÞ of the normalised Gibbs free-energŷ G are shown in Fig. 3a for selected values of the foundation stiffnessĵ with wðĜÞ denoting the density of states (i.e. the fraction of total number of morphological microstates that have a normalised free-energy in the rangeĜ toĜ þ dĜ). Two key features emerge: (i) probability of low free-energy states decreases with decreasing foundation stiffness and (ii) the probability density functions become more peaked with decreasing foundation stiffness. The normalised homeostatic temperature 1=f associated with these distributions is plotted in Fig. 3b wheref fjG S À ! 0 j. Consistent with the more uniform distributions pðĜÞ for the stiffer foundations, 1=f increases with increasingĵ. These results can be understood in terms of the competition between cytoskeletal and elastic energy discussed in Sect. 2.4 as follows.
With increasing cell area, the concentration of bound stress-fibres increases and therefore the concentration of the unbound proteins reduces. This increases the entropy of the stress-fibre proteins and reduces their contribution toĜ, i.e. contribution from the cytoskeleton becomes more negative as seen from the corresponding probability distribution ofF cyto in Fig. 3c. On the other hand, the elastic energy of the cell increases with increasing area and this gives rise to a minimum free-energy of the cell (Fig. 2a). We shall subsequently show that cells with larger areas exert higher tractions on the foundation. However, for stiff foundations, these tractions introduce small elastic energies in the foundation and consequently the minimum system free-energyĜ min for spread cells on stiff foundations is relatively low. By contrast, these same spread configurations introduce large elastic energies in compliant foundations with the consequence thatĜ min of the system with a soft foundation is higher than that for a stiff foundation. This implies that the system with a stiff foundation will explore free-energy configurations with a higher G so as to compensate and maintain the average freeenergy to be equal to G S . A wider distribution pðĜÞ with a high 1=f and a mode at lowerĜ then ensues for stiff foundations. We emphasize that the homeostatic temperature 1=f is much greater than the thermodynamic temperature T for the high stiffness foundations and thus the homeostatic ensemble permits larger fluctuations than those allowed by the conventional statistical ensembles. The physical origins of these high fluctuations are the exchange of high energy species such as glucose between the cell and the nutrient bath. This exchange causes large energy fluctuations that ultimately give rise to the large observed variability in experiments as we shall proceed to show.
To illustrate the multiplicity of morphological microstates with the same free-energy, some selected configurations of the cells on theĵ ¼ 22 and 2880 foundations are included in Fig. 4a, b, respectively. All these morphological microstates haveĜ values in the very close vicinity of the mode of the distribution (Fig. 3a). For each configuration, we have also included the distributions of the stress-fibre concentrations as parameterised by Fig. 3 Predictions of the a probability density functions pðĜÞ of the normalised Gibbs free-energy, b the normalised homeostatic temperature 1=f and c probability density functions pðF cyto Þ of the normalised cytoskeletal free-energy for selected values of the normalised foundation stiffnessĵ. In b we have indicated that cells in suspension have a zero homeostatic temperature It is clear that even for a given fixed free-energy, the cells can attain a large diversity of cell shapes, areas and distributions of cytoskeletal proteins. Thus, even though we have only shown a very small sample of highly probable states, these results are very much in line with the diversity of observations in experiments. Moreover, in line with observations [12,15], the images of the cells in Fig. 4a, b indicate that both the cell area and the level of actin polymerisation as parameterised byN b increase with increasing foundation stiffness. In particular, cells on stiff foundations have a high concentration of polymerised stress fibres at the cell periphery and in filopodia-like structures that set the cell polarity. We now proceed to quantify the dependence of such observables on the foundation stiffness.

Foundation stiffness affects traction forces and cell spreading
We focus on two key observables that are widely reported in micro-post experiments: (i) the cell area and (ii) the traction forces. With this in mind, we define a normalised cell area asÂ A=A 0 , where A and A 0 are the deformed and undeformed cell areas, respectively, while the normalised traction at location The normalised average traction then follows aŝ Probability density functions pðÂÞ and pT avg À Á are included in Fig. 5a, b, respectively, for a range of foundation stiffnesses. Similar to pðĜÞ, pðÂÞ and pT avg À Á become more peaked with decreasing foundation stiffness with the mode of the distribution simultaneously shifting to a lowerÂ andT avg . Thus, in line with experimental measurements [15] we predict that not only do the observed cell areas decrease with decreasing foundation stiffness but also the increasingly peaked distributions with decreasingĵ imply smaller standard errors in measurements. The overall reason for this is similar to that discussed above whereby cells on stiff foundations can spread more to lower their free-energy without introducing a large elastic penalty from the foundation and thus can sample a wider variety of morphological microstates. We note in passing that another observable typically reported in experiments is the cell aspect ratio as defined by the ratio of the major to minor axes of a best fit ellipse. For the cases investigated here the aspect ratio distribution was relatively insensitive to the foundation stiffness with the cell having a mean aspect ratio % 2:7 in all cases.
The Markov Chain Monte Carlo calculations used in the simulations are able to determine the entire probability distribution functions (Fig. 5) by sampling in excess of few million equilibrium morphological microstates. On the other hand, experiments typically report statistics based on observations of 10-50 cell configurations [12][13][14][15] and are therefore unable to generate distributions of the type in Fig. 5. Rather, experimentalists commonly plot statistics in the form of so-called box-and-whisker diagrams. Hence to make more definitive contact with measurements, the data in Fig. 5 is re-plotted in Fig. 6a, b in the form of box-and-whisker diagrams for the distributions ofÂ andT avg , respectively. The box plots depict the median and first and third quartiles of the distributions of these observables calculated from the sample list generated by the MCMC while the whiskers show the 5th and 95th percentiles of the distributions. Clearly, the median values ofÂ andT avg increase with increasingĵ with the asymptote of the infinitely stiff foundation being approximately attained for j [ 3000. Moreover, the box-plots also clearly show that the diversity of observations increases with increasingĵ and this is most clearly seen in the wider spread of the whiskers with increasingĵ.

Correlation of traction forces and cell area
Consistent with measurements, the model predicts that the average tractionsT avg increase with increasing foundation stiffness. However, this rise inT avg is accompanied by enhanced cell spreading (i.e. larger cell area) with increasingĵ. Thus, it is unclear whether cells of a given area exert higher tractions on foundations with higher stiffness or if the increase seen in Fig. 6b is solely due to the increased cell area affectingT avg . Here we attempt to decouple these two effects.
Spatial distributions ofT for selected morphological microstates chosen from the mode of theĜ distributions are shown in Fig. 7a, b for foundations with stiffnessĵ ¼ 22 and 2880, respectively (these morphological microstates are the same as in Fig. 4). Consistent with the data in Figs. 5b and 6b and a host of measurements [12][13][14][15], we see that the cells exert smaller tractions on the more complaint foundations. Moreover, as clearly seen in Fig. 7b the tractions are generally highest along the cell periphery. Of course, cells assume a large diversity of shapes and sizes even on the mode of theĜ distribution (the diversity being  larger for the cell on the stiffer foundation) and therefore it is not possible to decouple the effects of cell size and foundation stiffness directly from Fig. 7. We thus examine the homeostatic ensemble as sampled by the Markov chain as follows. For a given cell area, the Markov chain has a large number of morphological microstates and thus we can determine the distribution of observables for a given cell area.
Here, we consider two observables, viz. the average tractionT avg and the nominal tractionT T T avg A=A 0 . The average and nominal tractions are equivalent to the average and total forces that are typically reported in experiments [12,15] of cells on arrays of microposts.
Predictions of the variation of the three quartiles (25th percentile, median and 75th percentile) ofT avg andT T are plotted in Fig. 8a, b, respectively, as a function of the spread cell areaÂ for cells on foundations with stiffnessĵ ¼ 22, 180 and 2880. For a given foundation stiffnessĵ, the dependence ofT avg onÂ is relatively mild (Fig. 8a). As a corollary, the nominal tractionsT T increase nearly linearly withÂ (Fig. 8b). However, it is clear from Fig. 8a that for a given cell areaÂ,T avg increases with increasingĵ and consequentlyT T too is higher for cells with the same area on stiffer foundations (Fig. 8b). We thus conclude that the increased average tractions exerted by the cell on stiffer foundations are primarily due to the formation of a more concentrated stress-fibre cytoskeleton in response to stiffer foundations, but is not due to larger cell size per se. Of course the stiffer foundations also cause increased cell spreading but while this results in the cell exerting a larger total force (i.e. nominal tractions), the increased cell spreading does not significantly increase the average tractions. Finally, we note that consistent with experiments [15], even for a given cell area there is a greater diversity in the tractions exerted by cells on stiffer foundations compared to the equivalent diversity for cells on more compliant foundations (note the logarithmic axis scale for tractions in Fig. 8). Thus, the larger diversity of tractions on stiffer foundations (Figs. 5b, 6b) is not solely due to the larger diversity of cell areas that the cell assumes on stiffer foundations. This is understood  by recalling that a morphological microstate is not solely characterised by a single observable such as cell area: the higher homeostatic temperature on stiffer foundations (Fig. 3b) implies that the cell assumes a great diversity of morphological microstates all with the same area and these microstates can exert a large range of tractions.

Concluding remarks
We have used the homeostatic mechanics framework of Shishvan et al. [22] to analyse the response of cells on a dense array of micro-posts idealised as a Winkler foundation. The framework enables the quantitative prediction of the stochastic response of cells with contractility, cell spreading and the traction forces that the cells exert on the foundation all coupled to the foundation stiffness. We show that in line with observations, the model accurately predicts that: (i) the extent of cell spreading, stress-fibre polymerisation and traction forces that the cells exert on the foundation increase with increasing foundation stiffness; (ii) the traction forces are primarily concentrated along the cell periphery; and (iii) while the total tractions increase with increasing cell area, the average tractions are reasonably independent of cell area, i.e. for a given substrate stiffness, the average tractions that are normalized by cell area do not vary strongly with cell size.
These results thus clarify the question on the source of the increased tractions with increasing foundation stiffness. In particular, they suggest that the stronger tractions that cells exert on stiffer foundations are a result of the higher levels of stress-fibre polymerisation and not directly linked to the larger levels of cell spreading that also occur on stiffer foundations.
A key feature of the model is that the predictions are statistical with cell shape, size, tractions and all other observables being outcomes of the predictions in the form of probability distributions. In fact, in line with observations, the framework predicts that the diversity of most observables such as cell area and tractions increases with increasing foundation stiffness. The homeostatic mechanics framework thus suggests that the variability in experimental observations is inherent to the homeostatic equilibrium of cells rather than being a result of in vitro experimental errors.

Conflict of interest The authors have no conflict of interest to declare
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.