Contact guidance as a consequence of coupled morphological evolution and motility of adherent cells

Adherent cells seeded on substrates spread and evolve their morphology while simultaneously displaying motility. Phenomena such as contact guidance, viz. the alignment of cells on patterned substrates, are strongly linked to the coupling of morphological evolution with motility. Here, we employ a recently developed statistical thermodynamics framework for modelling the non-thermal fluctuating response of cells to probe this coupling. This thermodynamic framework is first extended via a Langevin style model to predict temporal responses of cells to unpatterned and patterned substrates. The Langevin model is then shown to not only predict the different experimentally observed temporal scales for morphological observables such as cell area and elongation but also the interplay of morphology with motility that ultimately leads to contact guidance. Supplementary Information The online version contains supplementary material available at 10.1007/s10237-022-01570-9.


Introduction
Living systems are characterised by a wide variety of timescales governing key biological processes at different length scales. For example, on a macroscopic scale, it is well known that animal species have their own circadian rhythm (Aschoff 1979;Ohata et al. 1999) that governs their daily behaviours, such as feeding and resting. On a microscopic scale, cells follow a cell cycle, which is a sequence of cellular phases such growth and division (Johnson and Walker 1999; Kastan and Burtek 2004). These timescales are critical for cell development and are controlled by complex molecular pathways (Johnson and Walker 1999; Kastan and Burtek 2004;Stanewsky et al. 1998;Kume et al. 1999). Another example pertains to the in vitro spreading of adherent single cells on substrates. Experiments have demonstrated that the spreading and elongation rates differ in time by an order of magnitude (Kesavan et al. 2014;Nisenholz et al. 2014). The rates of spreading and elongation were also observed to be non-constant, which suggested that distinct temporal phases may exist in a cell's morphological evolution.
Another important timescale pertaining to in vitro experiments of adherent cells is associated with cell motility. Motility occurs via large and co-ordinated morphological changes in the cell, e.g. the treadmilling mechanism in adherent cells, whereby cells move by tugging on the substrate by exploiting focal adhesions (Balaban et al. 2001;Bell et al. 1984). Therefore, the motility and morphological change timescales are expected to be interlinked. It is well known that substrate properties such as stiffness (Hadden et al. 2017;Ulrich et al. 2009), chemical composition (Engler et al. 2004;Carter 1965Carter , 1967 and topology (Buskermolen et al. 2020;Chang et al. 2013;Ray et al. 2017) strongly affect single-cell morphological and motile behaviour. Experiments have shown that changing the substrate stiffness alters both the spreading timescales of adherent cells and the speed at which they explore the environment (Asano et al. 2017;Ghibaudo et al. 2008;Pelham and Wang 1997). Guidance provided by substrate anisotropy, i.e. contact guidance (Guido and Tranquillo 1993;Dunn and Heath 1976;Weiss 1945), also influences cell spreading and motile behaviour. In fact, cells elongate more and faster on adhesive channels (Huang and Donald 2015) and this is also accompanied by an enhanced alignment of the biochemical force generating machinery within the cell, viz. stress-fibres (Buskermolen et al. 2019). Measurements have also reported a faster directional exploration speed in confined settings (Pathak and Kumar 2012). However, physical insights into the interplay between morphological changes and motility for cells on anisotropic substrates are lacking despite their importance for understanding cell guidance.
The aim of this theoretical work is to understand the timescales governing phenomena, such as contact guidance, that emerge due to the coupling of morphological evolution with cell motility. We propose a novel framework that we call the Homeostatic Langevin Equation (HLE) that is an extension to the Homeostatic ensemble (Shishvan et al. 2018) developed to understand the long timescale or stationary responses of cells. The HLE recognises the non-thermal fluctuating response of the cells and is used to simulate the temporal evolution of isolated cells seeded on unpatterned and patterned substrates. A schematic of the simulation set-up for a cell seeded on a substrate with a fibronectin stripe of width W is shown in Fig. 1a along with a representative prediction in Fig. 1b. The cell seeded from suspension is circular, and the HLE will predict its coupled morphological evolution along with its motility as shown in Fig. 1b. Key experimental observables such as cell area, cell shape, distribution of cytoskeletal proteins such as stress-fibre distributions, shape of nucleus and focal adhesion distributions are outcomes of the simulations. These predictions are used to extract the timescales of morphological evolution and cell motility and thereby develop an understanding of the mechanisms and processes that result in the guidance of cells on patterned substrates. Fig. 1 Simulation set-up. a Single cell placed on a fibronectin stripe of fixed width W . The cell exchanges high-energy nutrients with the nutrient bath. The components of the cell modelled explicitly include an elastic nucleus and cytoplasm as well as the contractile stressfibres in their polymerised state in equilibrium with their unbound components within the cytoplasm. b Temporal evolution of an adherent cell on a fibronectin stripe. In the simulation, cells are modelled as two-dimensional bodies in the x 1 − x 2 plane lying on a substrate. All cells are seeded from their state in suspension and then spread and change shape while simultaneously exploring the substrate. The images in (b) are representative images from simulations with the scalebar of 2R 0 the diameter of the circular cell in its elastic resting state

Modelling formulation
The interaction of cells with environmental cues (e.g. mechanical, chemical) has a profound influence on cell response. For example, not only do morphological observables such as cell area and aspect ratio increase with increasing substrate stiffness (Kesavan et al. 2014;Nisenholz et al. 2014) but cell motility also exhibits a biphasic behaviour on stiffer substrates (Peyton and Putnam 2005;Pathak and Kumar 2011). While different modelling frameworks are used to understand the influence of environmental cues on cell morphology (Vigliotti et al. 2015;Shenoy et al. 2016;Bengassser et al. 2013) and cell motility (Stokes et al. 1991;Klank et al. 2017), the interplay of the underlying mechanisms remains elusive. Our aim here is to develop a unified framework with the objective of better elucidating the underlying mechanisms. We do this by extending the recently developed homeostatic mechanics framework (Shishvan et al. 2018), which predicts the distribution of morphological states that an adherent cell assumes during the interphase period of the cell cycle.
The mentioned framework lacks temporal information, and here, we extend the formulation to enable predictions of the temporal evolution of morphological observables alongside cell motility as a function of environmental cues.

A brief overview of the homeostatic mechanics framework
The homeostatic mechanics framework recognises that a cell is an open system which exchanges nutrients with the surrounding nutrient bath (Fig. 1a); see Shishvan et al. (2018) for further details. These high-energy nutrient exchanges fuel large fluctuations (much larger than thermal fluctuations) in the cell response associated with various intracellular biochemical processes. The cell uses these biochemical processes to maintain itself in the homeostatic state. Specifically, homeostasis is the ability of a living cell to remain out of thermodynamic equilibrium by maintaining its various molecular species at a specific average number that is independent of the environment (Gordon 2013). This average number is sustained over all the non-thermal fluctuations of the cell (at-least over the interphase period of the cell cycle and in the absence of any imposed shock such as starving the cell of nutrients). The implication is that over the fluctuations of the cell from any reference state, ⟨ΔN i ⟩ = 0 where ΔN i is the change in the number of molecules of species i from its reference value with ⟨x⟩ denoting the average of x over the ensemble of states sampled over the non-thermal fluctuations. These fluctuations alter the cell morphology, and each morphological microstate (cell shape, protein distribution etc.) has an equilibrium Gibbs free energy G = ∑ i i N i , where i is the chemical potential of species i . Using the Gibbs-Duhem relation, we then rewrite this in terms of the reference state as where now 0 i is the chemical potential of species i in the reference state and G ref is the equilibrium Gibbs free energy of the cell in its reference state. Upon employing the homeostatic constraint that ⟨ΔN i ⟩ = 0 , we have ⟨G⟩ = G ref , i.e. irrespective of the environment, the ensemble average Gibbs free energy is constant. This is a universal constraint that quantifies the fact that living cells maintain themselves out of thermodynamic equilibrium but yet attain a stationary state which is the homeostatic state. While the above constraint specifies the average state of the cell, it remains to determine the distribution of states the cell assumes that satisfies this average. As mentioned, cells do not remain in a single microstate but fluctuate between different configurations after seeding. Cells explore these different morphological microstates thanks to several biochemical processes, such as actin polymerisation and treadmilling. The objective of the homeostatic mechanics framework is to capture the different configurations by adopting the ansatz that the observed distribution of cell shapes is the one with the overwhelming number of microstates, i.e. the distribution that maximises the morphological entropy subject to the homeostatic constraint, i.e. ⟨G⟩ = G ref = G S , with G S being the free energy of a cell in suspension, and to any other geometrical constraints, such as confinement, imposed by adhesive patterns on substrates. The distribution of states is denoted as an equilibrium distribution, where "equilibrium" means in a homeostatic state.
In this study, cells are modelled as two-dimensional bodies in the x 1 − x 2 plane lying on a substrate (Fig. 2). In this two-dimensional context, the cell morphology is defined by the positional vectors (q) r i (i = 1, 2) of the q = 1, … , M control points used to specify the cell shape. Denoting the positions of the control points in morphological microstate (j) as Shishvan et al. (2018) demonstrated that the cell in homeostasis assumes a given morphological microstate (j) with probability where Z ≡ ∑ j exp(− G (j) ) is the partition function of the morphological microstates, and the distribution parameter emerges from the homeostatic constraint = G S . Thus, 1/ in (1) is referred to as the homeostatic temperature, and it sets the equilibrium distribution of morphological microstates of the cell (also referred to as the homeostatic ensemble) as an analogous quantity to the thermodynamic temperature of the canonical ensemble (Gibbs 1902). The 2D cell in microstate (j) is described by the q = 1, … , M positional vectors (q) r i (i = 1, 2) . Here, we describe the numerical scheme used to define these vectors. In the 2D context of cells on substrates, describing a morphological microstate reduces to specifying the position of all material points where the cell is in contact with the substrate. Thus, given the location of material points X i on the elastically undeformed configuration, we impose a displacement field u (j) i (X i ) to obtain the displaced coordinates of these points x (j) i = X i + u (j) i which defines the morphological microstate (j) . The definition of the morphological microstate is then completed in the 2D setting by assuming that these points on the cell are connected to material points on the substrate in the same location x (j) i . The cell is modelled as a continuum, and thus, u (j) i is a continuous field. We define u (j) i via Non-Uniform Rational B-Splines (NURBS) such that the morphological microstate is defined by M control points of coordinates (q) r 1 , (q) r 2 . In Fig. 2 Simplification of an in vitro cell to a 2-dimensional body lying on the x 1 − x 2 plane all the numerical results presented here, we employ M = 16 with 4 × 4 control points (q) r 1 and (q) r 2 governing the displacements in the x 1 and x 2 directions, respectively. We construct the NURBS by using third-degree polynomials in both the x 1 and x 2 directions and by choosing two knots both with multiplicity four located at the extrema of the interval. The choice of these parameters enforces restriction on the sampled morphological microstates, i.e. features with wavelengths below the characteristic wavelength of the NURBS are not captured and thus typically this wavelength is chosen to represent the minimum width of a filopodium expected for the selected cell type (Shishvan et al. 2018;Buskermolen et al. 2019). Given the displacement fieldu (j) i , we evaluate the Gibbs free energy G (j) using model described below. The cell is discretised by constant strain triangles of sizee ≈ R 0 ∕10 , where R 0 is the radius of the cell in its elastic resting state.
The numerical methods used for sampling via Markov chain Monte Carlo (MCMC) the morphological microstates and to obtain 1/ are the same as Shishvan et al. (2018) and Buskermolen et al. (2019), and a brief description is provided in Appendix A.

Gibbs free energy of a morphological microstate
The implementation of the homeostatic mechanics approach described above requires a specific model for the Gibbs free energy of the cell-substrate system in a given morphological state. Modelling all the elements of the cell is unrealistic, given that many of their kinetics are still unknown, and often not required, as specific components are known to determine and control the cell response to different environmental cues. Here, we are interested in investigating the response of cells to adhesive patterning of the substrates. These cues are known to guide single-cell behaviour resulting in significant cell alignment bias as well as remodelling of the stress-fibre cytoskeleton. Thus, we use a model that includes stress-fibre cytoskeleton introduced by Vigliotti et al. (2016) and subsequently modified in Shishvan et al. (2018) as well as passive elasticity of the cytoplasm and nucleus. The Gibbs free energy of a morphological microstate (j) of the system comprises contributions from the cell and the substrate and is written as For the sake of brevity, we drop the superscript (j) throughout this section with the understanding that we are always discussing a given morphological microstate. The system (i.e. cell plus substrate) is within a nutrient bath and under atmospheric pressure as illustrated in Fig. 2. Then, taking atmospheric pressure as the reference zero state, G is given by as T e i = 0 on the surface S T of the system exposed to the environment with F cell and F sub the corresponding Helmholtz free energies of the cell and substrate, respectively. For the linear substrate model, the substrate Helmholtz free energy F sub = −G sub . The free energy of the cell comprises four contributions such that where f cyto , f ec , f passive are the cytoskeletal, energy carrier and passive Helmholtz free energies per unit volume in the cell, while f adh is the adhesive protein Helmholtz free energy per unit area over the interface S I between the cell and substrate. We define and the corresponding Gibbs free energy of the active components is given by the Legendre transform as where the active stress A ij is related to the strain ij via A ij ≡ (f cyto + f ec )∕ ij , while the surface tractions T i and stretch Δ i of the adhesion proteins are related by T i ≡ f adh ∕ Δ i . The free energy associated with the energy carriers drives mechanical work and so we assume and substituting (7) into (6) gives The Gibbs free energy of the system is then as follows: We assume that the collagen coating of both the unpatterned and patterned substrates is dense, and thus, there is not any impedance in the formation of focal adhesions; the term ∫ S I g adh dS is assumed constant for ever configuration. (3) (9) Additionally, the underlying substrate is assumed to be rigid, i.e. F sub = 0 . We proceed to discuss the formulation for g cyto and f passive .

Calculation of cell free energy
Calculation of g cyto and f passive involves modelling three elements within the cell: (i) a passive elastic contribution from elements such as the cell membrane, intermediate filaments, and microtubules in the cytoplasm, (ii) an active contribution from the contractile actomyosin stress-fibres that are modelled explicitly and (iii) the nucleus modelled as a passive elastic body. These are modelled using the framework of Vigliotti et al. (2016) and the subsequent modifications of Shishvan et al. (2018) and Buskermolen et al. (2019). We shall first describe the modelling of the active actomyosin stress-fibres in the cytoplasm and then discuss the elastic model of both the nucleus and the cytoplasm.
Consider an incompressible two-dimensional (2D) cell (both the cytoplasm and nucleus are assumed to be incompressible) of thickness b 0 , radius R 0 and volume V 0 in its elastic resting state comprising a nucleus of volume V N and cytoplasm of volume V C such that V 0 = V N + V C (Fig. 1a). We assume a cylindrical representative element (RVE) of volume V R = b 0 n R 0 ∕2 2 , where 0 is the length of a stress-fibre functional unit in its ground state and n R is the number of these ground-state functional units within the cytoplasm. The average number of functional units (polymerised and unbound) in the RVE is then N 0 = N T 0 V R ∕V C , where N T 0 is the total number of functional units present in the cell. The state of the stress-fibres in the RVE located in x i is described by their angular concentration x i , and the number of functional units n(x i , ) in series along the length of each stress-fibre. The angle is measured in the undeformed configuration and measures the orientation of the stress-fibre bundle with respect to the x 2 − direction of the stripe (Fig. 1a). At steady state, Vigliotti et al. (2016) showed that the number n ss of functional units within the stress-fibres is where ∼ ss nom is the strain at steady state within a functional unit of the stress-fibres, and (x i , ) is the cell stretch in direction , i.e. component of the left stretch tensor in direction . The chemical potential of the functional units within the stress-fibres is given in terms of the Boltzmann constant k B and of the normalised number N L of lattice sites available to unbound proteins by where N u ≡ N u ∕N 0 is the normalised concentration of the unbound stress-fibre proteins in terms of the concentration N u of the unbound, ̂ ≡ n R ∕N 0 and ̂ max is the maximum normalised value of ̂ corresponding to full occupancy of all available sites for stress-fibres (in a specific direction).
Here, the enthalpy b of n R bound functional units at steady state is given in terms of the isometric stress-fibre stress max and the internal energy b0 as where Ω 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 u of the unbound proteins as For a fixed configuration of the 2D cell (i.e. a fixed stretch distribution (x i , ) ), the contribution to the specific Gibbs free energy of the cell from the stress-fibres then follows as where 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 g cyto as N u (x i ) and ̂ (x i , ) are unknown and their value will follow from the chemical equilibrium of the cell as we shall discuss subsequently.
The total stress Σ 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. The total Cauchy stress is written by an additive decomposition as where A ij and 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 H 0 of the stress-fibre proteins as where is the angle of the stress-fibre measured with respect to x 2 , in the current deformed configuration. The passive elasticity in the 2D setting is given by a 2D specialisation of the Ogden (1972) hyperelastic model as derived in (Shishvan et al. 2018). The strain energy density function of this 2D Ogden model is for the cytoplasm and for the nucleus where I and II are the principal stretches, C ( N ) and C ( N ) the shear modulus and inplane bulk modulus of cytoplasm (nucleus), respectively, while m C ( m N ) is a material constant governing the nonlinearity of the deviatoric elastic response of cytoplasm (nucleus). The term P is a penalty term introduced to include an elastic penalty that prevents significant compression of the cytoplasm/nucleus in the x 1 − x 2 plane and given by P ≡ H J c − I II I II − J c 2 , where H(•) is the Heaviside step function, J c a non-dimensional constant that sets when this term becomes nonzero and sets the magnitude of the penalty. The cell is assumed to be incompressible, and thus throughout the cell, we set the principal stretch in the x 3 − direction III = 1∕( I II ) . The passive Cauchy stress then follows as i in terms of the principal passive Cauchy stresses p k ( ≡ k Ψ C ∕ k for the cytoplasm and ≡ k Ψ N ∕ k for the nucleus) and the unit vectors p (k) j (k = I, II) denoting the principal directions . The passive Helmholtz free energy of the cell is then f passive = Ψ C in the cytoplasm and f passive = Ψ N in the nucleus.
The equilibrium of a morphological microstate reduces to two conditions (Shishvan et al 2018): (i) mechanical equilibrium with Σ ij,j = 0 throughout the system, and (ii) intracellular chemical equilibrium. This second condition implies that the chemical potentials of bound and unbound stress-fibre proteins are equal throughout the cell and N u , by combining (11) and (13), as and N u follows from the conservation of stress-fibre proteins throughout the cytoplasm, viz.
Thus, knowing N u and ̂ x i , , the stress Σ ij can now be evaluated and these stresses within the system (i.e. cell and substrate) need to satisfy mechanical equilibrium, i.e. Σ ij,j = 0 . In this 2D case, the mechanical equilibrium condition is readily satisfied as the stress field Σ ij within the cell is equilibrated by a traction field T i exerted by the substrate on the cell such that bΣ ij,j = −T i , where b(x i ) is the thickness of the cell in the current configuration. Then, F cell becomes Here, u is given by Eq. (13) with the equilibrium value of N u obtained from Eqs. (19-20).

Morphological and other observables
The cell is modelled as a two-dimensional body whose morphology and position evolves on the substrate (Fig. 2). Typical observables such as 2D projections of the cell morphology and protein distributions can readily be extracted from the simulations.

Morphological observables
The morphological observables describe the cell shape and orientation in the unpatterned or patterned substrate. These are equivalent to analysing the 2D projection of a single adherent cell on the substrate. The observables of interest, as shown in Fig. 3, are: cell area A , cell aspect ratio A S , orientation and cell form factor FF . Note that all these observables are measured as a function of time.
Cell spread area A is calculated as the enclosed area within the outer boundary of the cell. The cell aspect ratio A S and orientation are defined from drawing the best fit ellipse of a given cell configuration. Note that cells take configurations that are not ellipses, but a best fit ellipse provides a metric to approximately define cell elongation and orientation and is extensively used in experiments (Buskermolen et al. 2019; Kesavan e al. 2014; Nisenholz et al. 2014; Ghibaudo et al. 2009). The best fit ellipse is calculated as an ellipse that best fits the outline of the current cell configuration. The cell aspect ratio is given by the ratio of the major to minor axis of the best-fit ellipse, and the cell orientation is the angle the major axis makes with the stripe direction. The cell form factor FF ≡ p 2 ∕(4 A) , where p is the cell perimeter, is a metric typically used to describe the elongation and capability of the cell to form filopodia on a given substrate. FF depends on the area but gives completely different information: FF = 1 means that the cell is circular in shape while FF ≫ 1 for a starlike cell, forming many filopodia to sense and explore the environment. The perimeter p is calculated as the perimeter of the polygonal shape drawn by the nodes on the periphery of the current cell configuration.

Visualisation of protein distributions and determining corresponding metrics
We show immunofluorescence-like images comprising of the nucleus, actin and vinculin stained in blue, green, and magenta, respectively (Fig. 3). This scheme is used for cell configurations shown throughout the paper as well as in the Supplementary Videos 1, 2 and 3. The nucleus is simply determined from colouring the mesh elements labelled as "nucleus" in blue. We proceed to discuss the protocol used to generate the protein staining for actin and vinculin. Actin stainings are representations of the stress-fibre distributions explicitly modelled in the framework. These distributions are plotted to convey two crucial pieces of Fig. 3 Illustrations of the cell morphological observables and staining distributions from a given morphological microstate. Scalebar is 2R 0 information: (i) the concentration of bound stress-fibre units in each location within the cell configuration and (ii) the local orientation of the densest stress-fibre bundle. While the formulation is a continuum formulation this information about stress-fibres is extracted from internal state variables and presented as follows. The concentration of stress-fibres ̂ n ss at each location x i and at an angle with respect to the x 2 − axis is known for any cell configuration. We then sketch an actin fibre in the direction max of the maximum value of ̂ n ss at that location and at a spacing that scales inversely with to convey that the actin network is more dense in that location. Focal adhesions are not explicitly modelled in the simulations. However, it is well known that mature and larger focal adhesions apply larger tractions on the substrate (Werner et al. 2017;Pelham et al. 1997;Frey et al. 2006). Thus, the magnitudes of the traction distribution T x i are used as surrogates to visualise the localisation of the focal adhesions. Specifically, higher traction magnitudes are represented by a deeper magenta colour to visualise suggest adhesion protein density in the given location.

Spreading of cells on substrates
It is instructive to briefly review the predictions of the homeostatic mechanics framework prior to extending the framework to predicting temporal evolutions. Recall that the homeostatic framework predicts the stationary distribution of microstates that the cell attains in a given environment. A cell in suspension needs to self-equilibrate and therefore assumes a unique configuration. For the 2D cell modelled here with the choice of parameters representing fibroblasts (Table 1), the cell in suspension is a circle of radius 0.92R 0 with the nucleus remaining undeformed with a radius R N . In this configuration, the elastic stresses generated by the compression of the cytoplasm are balanced by the tensile stresses generated by a spatially uniform distribution of stress-fibres within the cytoplasm.
Now consider a cell seeded on a patterned substrate as shown in Fig. 1b. The substrate is micropatterned with fibronectin stripes of a given width W , and adhesion of cell is prevented outside the stripes. This theoretical setup is equivalent to experimentally microprinting widely spaced adhesive stripes on the substrates such that the cell cannot span two adhesive stripes and is thus confined to a single stripe. The cell stresses no longer need to be self-equilibrated as intracellular stresses can be balanced via tractions between the cell and the substrate. The cell therefore no longer attains a unique configuration but rather can fluctuate between different morphological  Numerical penalty modulus J c 0.6 Critical value at which the elastic penalty starts to operate Shishvan et al. (2018) microstates with (1) specifying the probability of a specific morphological microstate and the distribution parameter set by the requirement that homeostatic constraint is satisfied over the entire ensemble of states the cell can attain. Predictions, obtained via MCMC (Appendix A), of the probability distribution of the normalised free energy Ĝ ≡ G∕|G S | of the cell are included in Fig. 4a along with the corresponding distributions of the cell area and aspect ratio in Fig. 4b and 4c, respectively, for three choices of the normalised stripe width Ŵ ≡ W∕(2R 0 ) . These metrics, typically used to characterise cell morphology, are defined as the normalised cell area A∕A R where A and A R are the areas of the cell and the area of the cell in suspension ( A R = 0.92R 0 2 ), respectively. Clearly, the cells are predicted to spread and elongate on the substrate with the driving force for these morphological changes arising from the fact that these spread and elongated states have low free energy: as discussed by Shishvan et al. (2018), spreading of the cell reduces its cytoskeletal free energy but results in an increase in the elastic energy and it is this competition that controls cell spreading. We observe that the confinement imposed by the adhesive stripes has a strong tendency to elongate the cells but has a minimal effect on the distribution of the cell free energy and cell area. It is important to recognise here that cells do not uniquely attain their lowest free-energy state but rather sample a wide distribution of morphological states as evidenced in Fig. 4b and 4c. This wide sampling is set by the distribution parameter as required to satisfy the homeostatic constraint. In fact, the model predicts that the normalised homeostatic "temperature" 1∕̂ = 1∕(G S ) increases with decreasing Ŵ as seen in Fig. 4d. These results are consistent with the predictions in Buskermolen et al. (2019) and serve as a useful reference as we proceed to understand the temporal evolution of cells seeded on substrates.

A Langevin style framework for cell dynamics
The homeostatic mechanics framework gives the stationary distribution (1) of morphological microstates that cells will attain in a given environment as seen above. This stationary state is typically attained within 24 to 36 h after seeding cells into that environment. The framework gives no temporal information about the evolution of cells as they attain this stationary state including whether all morphological observables reach their stationary distribution at the same rate or whether there are multiple operative timescales. Here, we propose the simplest possible extension of the homeostatic mechanics framework to a dynamical setting by invoking an analogy with the canonical ensemble. Within the context of the homeostatic ensemble, the Gibbs free energy G (j) of a cell fluctuates while the cell is in its stationary state (also referred to as homeostatic equilibrium), but the corresponding homeostatic potential M ≡ G − (1∕ )S Γ (Shishvan et al. 2018) remains constant, where G ≡ ⟨G (j) ⟩ and S Γ the morphological entropy of the cell. Thus, a direct analogy can be made between the homeostatic ensemble and the well-established canonical ensemble where the Helmholtz free energy F ≡ U − TS is a constant in a bath with a thermodynamic temperature T . In this bath, the system fluctuates over its microstates (j) such that it has a fluctuating internal energy U (j) which at equilibrium achieves average value U ≡ ⟨U (j) ⟩ and entropy S . Thus, we see that the internal energy U (j) and temperature T in the canonical ensemble are analogous to G (j) and 1∕ , respectively, in the homeostatic ensemble. The temporal evolution of the microstates of an isothermal system whose equilibrium is given by the canonical ensemble is often described by Langevin dynamics. The low speeds at which cells move (Purcell 1977) imply that it suffices to ignore inertia and overdamped Langevin dynamics for the homeostatic ensemble, which we call Homeostatic Langevin Equation (HLE), can be developed as follows. We specify that the temporal evolution of the co-ordinates r that describe the cell morphology is given by where t is time and is a damping coefficient, sometimes referred to as the mobility, that relates the velocity of the microvariable (q) r i to a determinstic force and f i (t) a contribution from a random force. Thus, while the first term on the right-hand-side of (23) represents a deterministic force contribution, the second term is the random contribution associated with large non-thermal fluctuations due to high-energy nutrient exchanges between the cell and the nutrient bath in which it resides. In principle, f i (t) are correlated in time since the molecular processes from which they originate have a finite correlation times (Roberts et al. 2014). However, observations (Stokes et al 1991;Roberts et al. 2014) suggest diffusive type behaviour of cells over timescales of a few minutes and thus we assume that there exists some correlation timescale decays rapidly. Then, on timescales ≫ t c (which in the case of fibroblasts is on the order of a few minutes), f i (t) can be taken to be a delta correlated stationary Gaussian process satisfying ⟨f i (t)⟩ = 0 and where ij and (• ) are Kronecker and Dirac deltas, respectively, and ξ the standard deviation of the random and decorrelated f i (t). Under these assumptions, (23) reduces to the Langevin equation. Thus, the HLE is valid for observation times scales ≫ t c (i.e. a vanishing Deborah number) and only captures the diffusive nature of cell motility but not the ballistic response of cells which occurs on the scale of a few minutes. It now remains to set the standard deviation ξ . In order to set ξ , we turn to the Fokker-Planck equation corresponding to the HLE. In the HLE, we expressed the uncertainty in the morphology of the cell in terms of Gaussian correlation functions. Shifting perspective, we can ask: what probability distribution P r, t,r 0 , where r 0 is r at time t = 0 , would give the same correlation functions? It is important here to stress that we do not care about that path r(t) the cell took but rather ask the simpler question of the probability that the cell attains a morphological microstate r(t) at time t , regardless of how it got there. This is equivalent to saying that if we run a large number of independent but nominally identical experiments then what is the probability distribution of morphological microstates observed at time t . This is a well-established problem and following Ichimaru (2018), the required probability distribution can be shown to satisfy the Fokker-Planck equation The steady-state solution to (25) corresponding to P(r, t)∕ t = 0 is the equilibrium probability distribution and given by: where Therefore, the Fokker-Planck Eq. (25) converges to the homoeostatic ensemble (1) by setting: We shall thus use this choice of ξ to resolve the temporal response of a single cell in a given environment. Recalling that f i (t) in (23) is random and decorrelated with a standard deviation ξ we can rewrite (23) in normalised form as: 1) is a Gaussian distribution of zero mean and unit variance. In writing (28), we have used (27) and the fact that the stochastic differential Eq. (28) is solved with a finite time step Δt withΔt ≡ Δt|G S |∕( R 2 0 ) . It is now apparent that while sets the fluctuation magnitude, the mobility sets the evolution timescale. We emphasise that this single timescale for the evolution of the morphological microstate does not imply that all observables evolve at the same rate as we shall proceed to show. All results are presented for model parameters representative of myofibroblasts as calibrated in Buskermolen et al. (2019), but we emphasise that the model holds more generally for any adherent single cell type.

Temporal integration of the HLE
The stationary distributions for the observables as well as the homeostatic temperature used in the dimensionless form of the HLE (28) were provided by MCMC calculations; see Appendix A. Stochastic differential equations are typically integrated using a Euler forward integration algorithm (Milstein 1994) which here was employed as follows: (i) The cell is placed from suspension on the origin of the substrate at time t = 0 . On the adhesive stripes the origin coincides with the axis of the stripe. (ii) For the microstate at time t described by the vectors (q)r i (t) , we evaluate the gradient Ĝ ∕ (q)r i using the function "Derivest" (D'Errico 2007). Derivest is a complex numerical differentiation algorithm that requires 46 separate evaluations of � G(r) ; readers are referred to D'Errico (2007) for details of the algorithm.
(iii) For all the degrees of freedom we then randomly select the value of the noise from the normal Gaussian distribution, N(0, 1) and insert it into (28) to obtain the value of (q)r i ∕ t . (iv) We then update the current configuration by via the forward Euler scheme.
(v) All the observables associated to the new cell configuration (q)r i t + Δt are then stored. (vi) Set t =t + Δt , and repeat from step (ii) until t =T Sim , where T Sim is the maximum simulated time of the experiment.
Time evolutions in both the unpatterned and patterned substrates were computed by selecting Δt = 0.001 . Numerical convergence checks confirmed that reducing the time step further did not result in any changes to the predictions. Moreover, at steady state all observables predicted by the HLE using this value of the time-step coincided with the MCMC calculations. All algorithms, except for Derivest, were developed inhouse in Matlab and the typical computational duration for a single HLE trajectory of duration T Sim = 4000 using 8 parallel CPUs is around 20 days.

Cell area and elongation evolve at different rates on unpatterned substrates
We shall discuss two cases: (i) cells on unpatterned substrates where there is effectively no confining effect imposed by the substrate within the x 1 − x 2 plane and then proceed to contrast with the case of (ii) contact guidance on substrates patterned with adhesive stripes of width W (Fig. 1b). We conceive of a typical experiment where a cell in suspension is seeded on a rigid substrate coated uniformly with an adhesive protein such as fibronectin. The evolution of the cell morphology is then observed as a function of time t , with t = 0 corresponding to the instant of seeding. Of course, motility of the cell and the evolution of its morphology are coupled, but first we focus on cell morphology. Representative images of computed cell morphologies and the corresponding actin, nucleus and focal adhesion organisations are included in Fig. 5a at three normalised times t = 1, 100, 3000 in addition to the state at t = 0 . Recall that the Langevin style dynamical Eq. (28) is a stochastic differential equation so that a different solution is generated for every realisation of the noise process, i.e. much like in experiments a different trajectory of morphological evolution is obtained for every solution of (28) with the same initial state at t = 0 . Hence, in Fig. 5a we show solutions at three times from three such trajectories. While the three morphologies computed using different trajectories of (28) are different, they show many similar features. These include the observations that with increasing time: (i) the cells spread and increase their area as well as their ellipticity or aspect ratio and (ii) the level of actin polymerisation and focal adhesions also increase.
It is instructive to make quantitative predictions of some morphological observables to enable more direct comparison with measurements. The complete morphology of the cell in the model is described by r (j) , but as in observations (Buskermolen et al. 2019) we focus on three coarse-grained metrics of the morphology, viz. (i) cell area A , (ii) cell aspect ratio A S and (iii) the form factor FF . Each solution of the HLE (28) produces a different trajectory, and much like in experiments, trends are most clearly seen by examining the ensemble average over a number of trajectories. We performed n = 100 independent simulations and defined the ensemble average as follows. At time t , the ensemble average area is defined as is the area of the cell in the i th trajectory at time t . The ensemble average aspect ratio A S (t) and FF(t) are defined in an analogous manner. While A S and FF are non-dimensional, it is instructive to define a nondimensional area as where A ∞ ≡ A(t → ∞) so that Â assumes values of 0 when the cell is seeded and fluctuates around 1 at convergence. Predictions of the temporal evolution with normalised time t Fig. 5b and are consistent with the qualitative features seen in images of the cell morphologies (Fig. 5a). However, an intriguing feature emerges: while the HLE has only a single timescale, the normalised cell area Â evolves and reaches its steady state much faster compared to the cell aspect ratio and form factor, i.e. different timescales for the evolution of area and shape emerge from the HLE. See Supplementary Video 1 for the simulated evolution of the cell over a single trajectory shown alongside the plots of the evolution of the ensembled averaged observables. The different timescales for the evolution of cell area and aspect ratio have also been observed in experiments. In Fig. 6a and 6b, we reproduce measurements from  Kesavan et al. (2014) of the temporal evolution of normalised area and aspect ratio, respectively, for the spreading of fibroblasts on stiff substrates with t = 0 corresponding to the instant of seeding the cell from suspension. The experimental data for area were normalised following (30) while the aspect ratio was inferred from the circularity data reported by Kesavan et al. (2014). In Fig. 6a, the average cell area over about 50 measurements is plotted with the error bars corresponding to the standard deviation, while in Fig. 6b the average aspect ratio is plotted over 10 measurements (two measured trajectories are also included to give an indication of the large variability in this metric over different trajectories). The striking difference between Fig. 6a and 6b is that while the steady-state cell area is achieved in about 150 min the cell aspect ratio evolves significantly more slowly such that is unclear whether a steady-state aspect ratio is achieved 20 h after seeding. We superimpose on Fig. 6 our corresponding predictions (averaged over 100 different trajectories) with the timescale of in the simulations chosen to be ( R 2 0 )∕|G S | = 10min so as to match the simulated and measured timescales of the temporal evolution of cell area. With this choice of the timescale we see that the simulations capture the observation that the cell aspect ratio evolves significantly more slowly. Moreover, the simulations not only capture the average temporal evolution response but also agree with measurements of the standard deviation of cell area. In addition, two selected simulated aspect ratio trajectories are included in Fig. 6b and illustrate that, in line with measurements, the simulations also display a very large variability of this metric over different trajectories.
The success of the simulations in capturing the fact that the two morphological observables evolve at different rates suggests that a single rate constant ( R 2 0 )∕|G S | suffices to set the evolution of these observables. We therefore use the model to interrogate the source of the two timescales that set the evolution of cell area and aspect ratio. For this, we consider a significantly simplified model where the cell is restricted to remain a spatially uniform ellipse with a fixed orientation (Fig. 7a). In this case, the stretches 1 and 2 of the principal axes of the ellipse completely define the cell morphology rather than the M positional vectors (q)r i . The deterministic evolution of the cell morphology in this simplified model is then given by an equation analogous to (28) with the noise term neglected, viz.
where i = 1, 2 . Predictions of the temporal evolution of Â and A S (where cell area A = 1 2 R 2 0 and A S = 1 ∕ 2 with 1 ≥ 2 ) are included in Fig. 7b for a cell seeded from suspension at t = 0 . Intriguingly, the qualitative feature that the cell area attains its steady state significantly faster than aspect ratio is retained in this very simplistic setting.
We observe three temporal regimes: (i) an initial regime I of rapid cell spreading where the cell remains circular with A S = 1 ; (ii) a subsequent regime II of cell elongation where cell area is constant but the aspect ratio increases and (iii) a final regime III where both cell area and aspect ratio increase although the changes in this regime are relatively minor. Regimes I and II set the two observed timescales for the evolution of cell area and aspect ratio, respectively. In this deterministic and simplistic setting where the cell morphology is only a function of ( 1 , 2 ) , we can understand this by examining the energy landscape Ĝ 1 , 2 shown in Fig. 7c. We include in Fig. 7c  The simulations use a timescale ( R 2 0 )∕|G S | = 10min . This single timescale is shown to capture the observation that the cell area evolves in minutes, while the cell aspect ratio evolves over a timescale of many hours space starting from its state in suspension along with isolines of A∕A R = 1 2 ∕(0.92) 2 and A S = 1 ∕ 2 . The trajectory set by Eq. (31) has two distinct branches: (i) an initial branch corresponding to regime I where the cell traverses a path of A S = 1 while A increases and then (ii) turns and traverses a path of constant A but with increasing A S . This trajectory is purely set by the topology of the Ĝ ( 1 , 2 ) landspace and the fact that (31) requires the cell morphology to evolve along a path with the steepest gradient in Ĝ ( 1 , 2 ) space. Thus, we argue that the two different timescales for the evolution of cell area and aspect ratio are purely a result of the freeenergy landscape of the cell. This landscape is set by the interplay between the elastic energy and cytoskeletal energy of the cell.

Diffusive motility of cells on unpatterned substrates
When cells are seeded on substrates coated with an adhesive protein, not only does their morphology evolve but this morphological evolution is coupled to their motility. The HLE (28) enables predictions of this coupled motility and morphological evolution of the cell. Predictions of Fig. 7 A reduced model of a cell as a spatially uniform ellipse to illustrate morphological evolution. a Sketch of the cell including its nucleus as a spatially uniform ellipse. Here, 1 and 2 are the stretches of the axes of the ellipse. b Predictions of the temporal evolution of the normalised ellipse area and aspect ratio A S ≡ 1 ∕ 2 ( 1 ≥ 2 ) for a cell seeded on an unpatterned substrate from suspension at time t = 0 . The predictions are for a deterministic response with effects of biological noise excluded. c Free-energy landscape of the reduced cell model. Contours of the normalised free energy Ĝ are included on map with axes 1 and 2 along with contours of A∕A R and A S . On the map, we show the deterministic trajectory of the cell in the freeenergy landscape for a cell starting from its state in suspension until it attains its minimum free-energy state on the substrate the coupled evolution of the cell morphology and its motility are shown in Fig. 8 for 4 selected trajectories of (28) with the cell seeded from suspension at same location at time t = 0 in each case. Stochastic motility (i.e. cells take a random path over the substrate surface) is predicted over the timescales in Fig. 8a in line with numerous observations (Stokes et al. 1991;Plou et al. 2018;Dunn and Brown 1987;Schienbein and Gruler 1993;Krummel et al. 2016). The stochastic motility can be characterised in terms of the squared displacement of the centroid of cell (k) from its seeding location (taken to be the origin of the Cartesian co-ordinate system with x (k) (t) and y (k) (t) the Cartersian co-ordinates of the cell centroid at time t , where x ≡ x 1 and y ≡ x 2 ). Predictions of R 2 (k) ≡ R 2 (k) (t)∕R 2 0 are included in Fig. 8b for 5 trajectories: we predict large variations in R 2 (k) over different trajectories consistent with observations (Stokes et al. 1991). (The timelapse of the motility of a single trajectory of a cell along with the corresponding R (t) is included in Supplementary Video 2.) Given this stochastic nature of the motility, it is more instructive to consider the mean squared displacement (k) over n independent trajectories. The predictions of M SD ( n = 100 ) are included in Fig. 8b and illustrate the diffusive nature of the predicted motility as M SD scales linearly with normalised time t . If the cell is assumed to be a single particle whose position is given by the cell centroid, the HLE (28) would predict a diffusion coefficient D = 1∕( ) such that the slope of the M SD curve in Fig. 8b is given by the normalised diffusion coefficient D ≡ D ∕|G S | = 1∕( |G S |) . (The HLE in this context describes the Brownian motion of the cell with the homeostatic temperature setting the magnitude of the fluctuations.) The M SD slope in Fig. 8b ≈ 0.06 , while for the cell parameters used here D ≈ 0.23 , i.e. the actual motility of the cell as predicted by (28) is significantly more sluggish compared to modelling the cell as a Brownian particle with a diffusion coefficient D . The slower motion is related to the morphological evolution of the cell that is coupled to its motility. If the cell were modelled as a single Brownian particle, all fluctuations would result in motion of the particle. On the other hand, in (28), the fluctuations affect the positional vectors (q) r i (i = 1, 2) of the q = 1, … , M that describe the cell morphology and a displacement of the cell centroid requires the co-ordinated motion of these vectors. Such co-ordinated motion via stochastic fluctuations occurs with a low Finally, we emphasise that our model captures the coupled evolution of cell morphology and motility over timescales ≫ t c , i.e. we only capture the diffusive motion of the cell and not its ballistic motion (Dunn and Brown 1987) due to persistence in the direction of cell motion over timescales ≪ t c . A consequence of the model only capturing the diffusive part of the motion is that it does not provide insights into the physical processes, such as actin treadmilling, that give rise to ballistic motion (Roberts et al. 2014;Lee et al. 2021;Recho et al. 2013;Cardamone et al. 2011).

Fig. 9
Coupled motility and morphological evolution of the cell seeded on a substrate with adhesive stripes. a Two Langevin trajectories for a cell seeded at time t = 0 from suspension on a a substrate with adhesive stripes of normalised width Ŵ ≡ W∕(2R 0 ) . One trajectory shows the evolution up to time t = 250 while the second to t = 1500 . Results are shown for three stripe widths Ŵ = ∞, 4 and 1. The scalebar 2R 0 is the diameter of the circular cell in its elastic rest state b, c. The corresponding temporal evolution of (b) the normalised cell area Â and aspect ratio A S averaged over n ≥ 50 Langevin trajectories. d Dependence of the timescales of the evolution of the area and aspect ratio on Ŵ . Cells no longer remain attached to the stripes for � W < 0.85 (Buskermolen et al. 2019), and this region is shown shaded 1 3

Contact guidance on substrates with adhesive stripes
In vivo, adherent cells interact with the surrounding extracellular matrix (ECM) that is not only responsible for the structural integrity of tissues but also establishes and maintains the cellular microenvironment by providing cells with mechanical, biochemical and physical cues. It is now well established that cellular microenvironments induce cells to align and migrate along the direction of the anisotropy-a phenomenon called contact guidance (Chang et al. 2013). Various in vitro chemical micropatterning approaches using two-dimensional (2D) substrates have been developed to study cellular contact guidance, as model systems to simplify the highly complex three-dimensional (3D) environments in vivo (Buskermolen et al. 2019). A common micropattern is fibronectin stripes of a given width W with adhesion of cell prevented outside the stripes (Buskermolen et al. 2019;Buskermolen et al. 2020). We now proceed to investigate via the HLE the effect of the confinement provided by the stripes on the evolution of cell morphology and the accompanying motility which finally leads to contact guidance. Given a cell seeded in the middle of the stripes at time t = 0 from suspension, Fig. 9a shows snapshots of the cell morphology and position from a single trajectory on stripes of normalised width Ŵ ≡ W∕(2R 0 ) = 4 and 1 as well as the unpatterned substrate with Ŵ = ∞ . The snapshots are shown for two selected times along with best fit ellipses and the corresponding orientations of the cell with respect to the stripe. From the n ≥ 50 trajectories computed here we have selected in Fig. 9a trajectories where the cell motion backed onto itself between times t = 250 and 1500 so that the cell positions are well separated, and we can easily illustrate the cell morphologies. While clear visual differences between the cell morphologies are observed for the cell on the different stripes, the most striking differences are in the cell trajectories where we clearly see that the cell is "guided" on the Ŵ = 1 stripe with the trajectory being nearly onedimensional (also see Supplementary Video 3).
A key outcome of the discussion in Sect. 4 is that on unpatterned substrates the area of the cell evolves much faster than its aspect ratio. Is this feature carried forward to when cells are constrained to within a single adhesive stripe? To investigate this, we observe from Fig. 9b and 9c that the key morphological variables of interest, viz. the area and aspect ratio, evolve in approximately an exponential manner. Thus, we fit a curve of the form to their temporal evolution averaged over n trajectories.
Here, x denotes the value of the observable of interest (i.e. either normalised area Â or aspect ratio A S ) with x ∞ denoting the value at normalised time t → ∞ while x R is the value of the observable for the cell in suspension at t = 0. While x ∞ and x R are directly available from trajectories as shown in Fig. 9b and 9c, the time constant x is obtained via a best fit of (33) to the Langevin predictions. The inferred time constants Â and A S for the evolution of the area and aspect ratio, respectively, are included in Fig. 9d as a function of the stripe width Ŵ . (Cells no longer remain attached to the stripes for � W < 0.85 (Buskermolen et al. 2019), and this region is shown shaded.) With decreasing Ŵ the aspect ratio evolves faster while the time for the area to attain its steady-state value increases. Consequently, Â and A S converge towards each other at small Ŵ eliminating the difference in the timescales for the area and aspect ratio evolution seen for the unpatterned substrates. For cells seeded on stripes, even constraining ourselves to spatially uniform elliptical cells results in a complex energy landscape as not all values of the ( 1 , 2 ) are permissible for all cell orientations . Thus, a simple picture like that presented in Fig. 7 (33) Fig. 10 a Evolution of the order parameter Θ for different stripe widths. b, c Motility of cells on substrates with adhesive stripes. Predictions of the temporal evolution of (b) normalised mean square dis-placement M SD and c dimensionality Λ of the motility on stripes of normalised width Ŵ ≡ W∕(2R 0 ) . The cell is seeded from suspension at time t = 0 on the stripes cannot be presented here to explain the convergence of the timescales for the evolution of area and aspect ratio with decreasing Ŵ . However, the only physics in the model that sets the timescales is the topology of the energy landscape. Hence, we conclude that the changes in the free-energy landscape due to the confining effect of the narrow stripes guide the area and aspect ratio of the cell to evolve hand-in-hand.

Timescale for the development of guidance and non-diffusive motility of cells on substrates with adhesive stripes
From the evolution of the morphological observables, we have discussed the dynamic effect of guidance, but we have yet to characterise its properties. This guidance can be characterised in two ways: (i) via the cell orientation that is related to cell morphology and (ii) via the direction of motion. Here, we first focus on cell orientation and in the following section on guidance which is more directly a consequence of motility.
To characterise the cell orientation related to changes in cell morphology with stripe width, we plot in Fig. 10a the order parameter Θ defined as where (i) is the orientation of the cell in the i th trajectory at time t (see Fig. 9a for the definition of (i) as the inclination of the major axis of the best fit ellipse to the stripe direction). The order parameter Θ is defined such that Θ = 0 if (i) is uniformly distributed over all n trajectories, while Θ = 1 if (i) takes a unique value. Cell alignment as parameterised through Θ also increases with decreasing Ŵ although the large increases in Θ occur at around Ŵ ≈ 1 when the aspect ratio increases. The increase in cell alignment at steady state with decreasing Ŵ has been reported in (Buskermolen et al. 2019). In that study, the authors argue that alignment results from the fact that cells near the edge of the stripes are necessarily aligned with the stripes: this boundary effect is of course more prominent for narrower stripes and also with cell elongation. Consequently, cells are more aligned for smaller Ŵ in line with the HLE predictions. Since alignment is a result of cells sensing the stripe edges, it follows that alignment is a consequence of cells wandering over the stripes including towards the edge of the stripes and hence an outcome of cell motility. Thus, not only does Θ increase with decreasing Ŵ but also the time required for the steadystate value of Θ to be attained decreases with decreasing Ŵ . This timescale for the evolution of Θ is of course set by the motility timescales of cells on the stripes.
To quantify the motility of cell on stripes, we include predictions of the variation of M SD ( n ≥ 50) with t on (34) stripes of different widths in Fig. 10b. Unlike the case of the unpatterned substrate ( Ŵ = ∞ ), the M SD on stripes of width similar to cell size (i.e. the Ŵ = 0.85 and 4 cases shown in Fig. 10b) does not seem to vary linearly with t . To understand the effect of the constraint of the stripes on motility, we introduce a measure of the mean square displacement in the direction of the stripes, viz.
where the seeding location at time t = 0 is assumed to be the origin of the co-ordinate system ( y ≡ x 2 − aligned with the stripe direction as shown in Fig. 9a). To parameterise the influence of the constraint of the stripes, we then define This parameter that characterises the dimensionality of the motility is defined such that Λ = 1 on an unpatterned substrate as M SD = 2MSD y , since movements in the x ≡ x 1 and y− directions are free and non-distinguishable, while Λ → 0 for � W ≪ 1 when the motion of the cell is completely constrained to be only in the y− direction (i.e. is one-dimensional) so that M SD =MSD y . Predictions of the temporal variation of Λ are included in Fig. 10c for selected stripe widths Ŵ . In all cases, Λ ≈ 1 at early stages of the cell motion, but, except for the unpatterned substrate ( Ŵ = ∞) , Λ subsequently reduces to 0. However, the time at which this transition from 2 to 1D motility occurs increases with increasing Ŵ . This can be understood by recognising that the x− displacement of the cell is constrained by the stripe width while the y− direction displacement is unconstrained and thus with increasing time M SD is dominated by the y− direction displacement and M SD →MSD y . The transition from a two-dimensional (2D) motion of the cell during the early stages of cell motion to 1D motility in the later stages (with the transition time being Ŵ dependent) induces the loss of the linear dependence of M SD with time seen in Fig. 10b for the finite stripe widths. Moreover, we would also anticipate that the attainment of the steady state of the order parameter Θ is governed by the time of transition from 2D to 1D motility. A comparison of Fig. 10a with 10c confirms that indeed Θ attains it steady-state value at approximately the time that Λ → 0 and thus the time to achieve steady-state alignment decreases with decreasing Ŵ . Finally, we note that the transition of Λ from 1 to 0 denotes the guidance of the cells by the stripes with the time taken for Λ to transition from 1 to 0 the time for contact guidance to be achieved on the patterned substrate.

The dynamics of biochemical alignment are strongly dependent on the morphological evolution
So far we have discussed the effect of confinement on cell motility and evolution of morphological observables. The model also explicitly includes the stress-fibre polymerisation that are strongly linked to cell morphology and here present predictions of the HLE for the evolution of the cytoskeletal stress-fibre arrangements in terms of a cytoskeletal order parameter that we now proceed to define. Given the spatial distribution of the direction max of the dominant stress-fibre bundle (set by the maximum value of ̂ n ss ) at any location x i , we then define a rotationally invariant spatial distribution of stress-fibre orientations as Predictions of the spatial distribution of ̂ (x i ) for cells seeded on stripes of different widths and at selected times t in a given HLE trajectory are shown in Fig. 11a (also see Supplementary Video 4). Homogeneity of colour in Fig. 11a is an indication of stress-fibre alignment and while the results in Fig. 11a suggest that stress-fibre alignment increases with time, the more dominant effect is seen with decreasing stripe width. To quantify this, we define a stress-fibre order parameter Θ cyto analogously to (34) with replaced by ̂ . The temporal evolution of Θ cyto for Ŵ = ∞, 4, 1 (averaged over n ≥ 50 HLE trajectories) is shown in Fig. 11b. While the predictions for Ŵ = ∞ and 4 are nearly identical, a rather interesting feature is observed for the Ŵ = 1 case with Θ cyto increasing rapidly and then decaying to a steady-state value (37) ≡ max − 1 V C � V C max dV. Fig. 11 a Snapshots of the distribution of stress-fibre orientation ̂ for cells seeded on stripes of different normalised widths Ŵ at selected times t . b Evolution of the cytoskeletal order Θ cyto as a func-tion of t with a zoom-in showing the early time evolution. c Evolution of the aspect A S in the early stages corresponding to the zoom-in in (b). The cell is seeded from suspension at time t = 0 on the stripes that is higher than that for Ŵ = ∞. To better understand this rapid increase in Θ cyto , we show in Fig. 11b a zoom-in for t ≤ 100 where we observe that the maximum value of Θ cyto is attained at t ≈ 50 . To examine the relevance of this time, we include in Fig. 11c the corresponding temporal evolution of the aspect ratio where we observe a biphasic evolution: a rapid increase for t ≤ 50 followed by a slower increase. We hypothesise that the rapid increase in aspect ratio is accompanied by or even requires high stress-fibre alignment and subsequently the alignment relaxes to increase the entropy of the stress-fibres and decrease cell free energy. These complex couplings are natural outcomes of the HLE.

Concluding discussion
The temporal response of isolated cells on unpatterned and patterned substrates has been investigated via a novel framework, labelled the homeostatic Langevin equation (HLE), that recognises the non-thermal fluctuations of cells. The morphological evolution is driven by gradients of the cell free energy, and we show that the HLE correctly predicts that the cell area or spreading evolves at a rate an order of magnitude faster than cell aspect ratio or elongation. These two very different timescales emerge as a consequence of the interplay between the stress-fibre polymerisation and cell passive elasticity. The framework enables the prediction of the coupled evolution of cell morphology along with cell motility. Over the timescales when the HLE is applicable, the simulated cell motility on an unpatterned substrate is Brownian, in line with numerous observations (Stokes et al. 1991;Plou et al. 2018;Dunn and Brown 1987), and emerges from coordinated morphological fluctuations. On substrates patterned with adhesive stripes cells again spread and elongate much like on unpatterned substrates. However, key differences now emerge related to the fact that the Brownian motion of the cells is now restricted outside of the parallel direction to the stripes. For example, for cells seeded on stripes the timescales for the evolution of area and aspect ratio evolution converge with decreasing width of the stripes. Furthermore, for narrow stripes, the evolution of the stress-fibre cytoskeletal arrangement shows a biphasic behaviour with an initial rapid rise associated with a rapid increase in the cell aspect ratio and then a decay to its steady-state value. The differences in the dynamics of the morphological observables carry to cell motility too. The HLE predicts that if sufficient time is given for the cells to explore the stripe widths, their 2D Brownian motion reduces to one-dimensional (1D) motion along the length of the stripes. This results in an apparent non-Brownian effect with the mean-square-displacement of the cell centroid no longer scaling linearly with time over the entire duration. A more important consequence of this switch to 1D motility is contact guidance or rather 1D motion of the cell along the stripes which manifests itself also in terms of the alignment of the cell orientation with the stripe direction. Thus, a key conclusion of the HLE framework is that the non-thermal cell fluctuations give the cell its ability to explore the stripe width and in turn, rather counterintuitively, result in its guidance and alignment with the anisotropy of its environment. else repeat the configuration (j − 1) in the sample list. Return to step (ii). (vi) Keep repeating this procedure until a converged distribution is obtained. Typical Markov Chains consist of N = 1 million samples. (vii) After running k different Markov Chains (to ensure that the correct sampling of the phase space), calculate the ensemble Gibbs free energy ⟨G (j) ⟩ = 1 N ∑ N j=1 G (j) , where N = kN is the total number of sampled configurations over the k Markov Chains. If ⟨G (j) ⟩ is within ±2% of G S , then this distribution is accepted, else the value of is modified and the process is repeated from (i).
A "burn-in" number of initial microstates from each Markov chain are discarded in the chain to eliminate the dependency of the sampled configurations from the initial random guess. The value of burn-in we used is 50, 000 . As is typical in such Markov chain Monte Carlo (MCMC) calculations, we attempted to achieve an acceptance rate of about 35% in the Metropolis criterion and adjusted Δ to ensure that we stayed within ±5% of this target acceptance rate. All algorithms, including the free energy evaluations, were run using in-house developed Matlab codes, and the duration to compute a typical Markov chain is between 1 and 2 h.