## Abstract

A model for the contractility of cells is presented that accounts for the dynamic reorganization of the cytoskeleton. It is motivated by three key biochemical processes: (i) an activation signal that triggers actin polymerization and myosin phosphorylation, (ii) the tension-dependent assembly of the actin and myosin into stress fibres, and (iii) the cross-bridge cycling between the actin and the myosin filaments that generates the tension. Simple relations are proposed to model these coupled phenomena and a constitutive law developed for the activation and response of a single stress fibre. This law is generalized to two- and three-dimensional cytoskeletal networks by employing a homogenization analysis and a finite strain continuum model is developed. The key features of the model are illustrated by considering: (i) a single stress fibre on a series of supports and (ii) a two-dimensional square cell on four supports. The model is shown to be capable of predicting a variety of key experimentally established characteristics including: (i) the decrease of the forces generated by the cell with increasing support compliance, (ii) the influence of cell shape and boundary conditions on the development of structural anisotropy, and (iii) the high concentration of the stress fibres both at the focal adhesions and at the sites of localized applied tension. Moreover, consistent with the experimental findings, the model predicts that multiple activation signals are more effective at developing stress fibres than a single prolonged signal.

## 1. Introduction

Most living cells sense, support or generate forces. For example, skeletal and heart muscle cells generate contractile forces on excitation, thereby performing many essential body functions. Endothelial cells can recognize the magnitude and the mode (steady or pulsating) of shear flows and respond accordingly, either maintaining a healthy endothelium or leading to vascular diseases. These and many other examples (see Bao & Suresh 2003) indicate that mechanical forces are central to the functioning of living cells. To develop models of the functional performance of a cell, it is not sufficient to just address its passive characteristics (except in a few simple cases; Hochmuth 1987). Instead, as discussed by Wang *et al*. (1993), force generation by the cytoskeleton must be included. The relevant forces have been measured using a succession of approaches. The first used continuous polymer substrates to examine the deformations (Harris *et al*. 1981). Later methods improved the resolution by (i) increasing the compliance of the substrates (Burton & Taylor 1997) and (ii) micro-patterning islands to facilitate measurements in areas as small as single focal adhesions (Balaban *et al*. 2001). More recently, Tan *et al*. (2003) introduced a method in which cells are laid on a bed of micro-needles: the independent deflections of each of these needles yield the forces exerted by the cell.

While significant experimental progress has been made to measure the forces generated by cells, interpretation of these experimental measurements poses a particular challenge. The reason is that a significant proportion of the forces are generated/supported by the cell cytoskeleton, which undergoes remodelling or reorganization in response to mechanical perturbations. This raises the fundamental paradox: how can we measure the intrinsic mechanical properties of living cells if they react to our measurement tools. Thus, unlike conventional engineering materials, *experiments alone provide few of the intrinsic properties of the cells*. The measurements can only be interpreted fully in the light of appropriate constitutive models. The development of a model for the contractility of the cytoskeleton is the main aim of this article.

Contractile bundles composed of actin and myosin proteins form and dissociate in non-muscle cells1 in response to external mechanical or chemical stimuli. Typically, these contractile bundles are found in three forms: (i) circumferential belts in epithelial cells, (ii) stress fibres along the ventral surfaces of cells, and (iii) contractile rings that form at the equator of a cell during cytokinesis. Here, we focus mainly on the stress fibres, though the general modelling framework developed here holds for other types of contractile bundles as well.

Tan *et al*. (2003) and Roure *et al*. (2005) have sought to correlate the forces exerted by mammalian cells on an array of micro-needles with the organization of visible stress fibres (figure 1). The forces are obtained from the deflections of the posts and the stress fibres revealed by subsequent application of an actin-staining procedure. It is apparent from the image in figure 1 (and numerous other images) that many of the force vectors are inclined to the axis of the visible fibre bundles. Indeed, some are almost normal to the bundles and, often, the largest vectors are present at locations where no visible stress fibres exist. The implication is that the forces are induced by fibres on a much finer scale, not revealed by staining procedures of this type. A corollary is that a contractility model capable of characterizing the forces should emerge from *continuum level* considerations, rather than from ensembles of discrete fibres. The primary intent of this article is to present a detailed description of a continuum model for the evolution of the cytoskeletal structure. Some capabilities of the model have been elucidated elsewhere (Deshpande *et al*. 2006), using it to simulate the response of a square cell attached by springs to four supports. These simulations highlighted the role of the spring compliance and of asymmetries in stiffness on the dynamic evolution of the stress fibres. The present article includes additional simulations which, together with those of Deshpande *et al*. (2006), demonstrate its consistency with key features found in experiments.

Previous attempts at developing models for the cytoskeletal network in stationary cells (i.e. neglecting cell spreading and motility) have taken the perspective that the cytoskeleton is an interlinked structure of passive filaments (Satcher & Dewey 1996; Storm *et al*. 2005). When included, cell contractility has been modelled by simply prescribing a thermal strain to either a cell regarded as (i) an isotropic elastic continuum (Nelson *et al*. 2005) or (ii) a discrete set of elastic filaments (Mohrdieck *et al*. 2005). Such models neglect the biochemistry of the active apparatus of the cell that generates, supports and responds to mechanical forces.

A generalized model for the contractility must be capable of characterizing the basic interactions between the forces, the assembly and dissolution of stress fibres and the compliance of the substrate. Moreover, once calibrated, it must explain such effects as the strong influence of substrate compliance on the forces, the dependence on cell size of the forces exerted at its periphery, as well the influence of cell shape and boundary conditions on the orientations of the stress fibres. It will be demonstrated that the model developed here is capable of addressing all of these effects. At this initial stage, incorporation of focal adhesions into the model is deliberately avoided. This exclusion allows a straightforward formulation capable of replicating several important features found in cells. The intent is to introduce the focal adhesions in later developments of the model.

The following logic has been used in the organization of this article.

The salient biochemical processes are summarized and brought into a framework that enables construction of the constitutive model.

The (three) coupled events that emerge from this assessment, expressed mathematically, become the fundamental concepts underlying the model. Initially, the concepts are implemented for a single fibre bundle to illustrate that the model has the elements needed to reproduce the following experimentally established characteristics: (

*a*) the force induced by the stress fibres decreases as the attachments become more compliant and (*b*) the average force per post increases with increase in the number of posts (Tan*et al*. 2003).The model is extended to two dimensions and used to simulate four features found in experiments: (

*a*) the cell develops strong structural anisotropy under uniaxial isometric loading while it remains isotropic under biaxial isometric loading, (*b*) as the stress fibres form, they concentrate adjacent to the attachments and orient preferentially, (*c*) multiple activations are more effective at rearranging the cytoskeleton into stress fibres, and (*d*) a high concentration of stress fibres form immediately adjacent to a site of localized applied tension.

## 2. Key biochemical processes and experimental observations

### (a) Biochemical processes

Here, we describe the key biochemical processes governing those aspects of the cytoskeleton dynamics that generate contractile forces in *non-muscle* cells. The main aim is to introduce the key processes that motivate the proposed bio-chemo-mechanical model, rather than review the extensive literature on focal adhesions, contractility and signalling (see Appendix A in the electronic supplementary material (ESM) for further details and schematics of the various processes).

In the *suspended or resting state*, the binding proteins or integrins are dispersed over the cell surface (and may be attached to some actin filaments). The short actin filaments in the cytoplasm are surrounded by a pool of actin monomers bound to profilin. Myosin II exists in the bent state in which the tail domain interacts with the motor head. *The formation of stress fibres* in the cell is triggered by an activation signal in the form of either a nervous impulse or an external signal. Several parallel intracellular pathways are involved. For example, adhesion to the extracellular matrix triggers a signalling pathway that induces the activity of profilin, cofilin and gelsolin. In turn, this activates PLC, which hydrolyses PIP_{2} and stimulates the release of Ca^{2+} from the endoplasmic reticulum. The influx of Ca^{2+} activates gelsolin, which cleaves the capped actin filaments into tiny fragments. The large numbers of free ends generated in this manner are rapidly elongated by the monomeric actin pool, forming many long filaments, some cross-linked with filamin and some bundled by α-actinin. Phosphorylation triggered by Ca^{2+} causes myosin II to preferentially assume its extended state. This promotes the assembly of myosin II into bipolar filaments that enter into the α-actinin-bound actin filament bundles, resulting in the formation of stress fibres. These fibres generate tension by cross-bridge cycling between the actin and the myosin filaments. When the tension is allowed to relax, the actin filaments are no longer held in place by the bipolar myosin filaments and the stress fibres disassemble.

### (b) Key experimental observations

Substantial evidence supports the idea that tension contributes to the formation of stress fibres. One set expresses the following role of suspended fibroblasts in the contraction of collagen gels. Free-floating gels contract over several days, by as much 90%, even though they lack stress fibres (Burridge & Chrzanowska-Wodnicka 1996). If the gels are anchored, isometric tension is generated and the fibroblasts develop prominent stress fibres (Mochitate

*et al*. 1991). Upon release of the tension, the attached gels rapidly contract, followed by stress-fibre disassembly (Mochitate*et al*. 1991; Grinnell 1994). Additionally, application of tension to cells in culture stimulates the formation of stress fibres (Franke*et al*. 1984). Moreover, when force is applied locally, an actin filament bundle is induced immediately adjacent to its application site (Kolega 1986). Associated with the tension-dependent stress-fibre assembly is the development of structural anisotropy. For example, uniaxially constrained, fibroblast-populated collagen gels develop a high degree of fibre alignment and mechanical anisotropy, while gels constrained biaxially remain isotropic (Thomopoulos*et al*. 2005).Cells precisely sense the restraining force and respond by localized proportional strengthening of the cytoskeleton linkages, allowing a stronger force to be exerted on the integrins (e.g. Choquet

*et al*. 1997). This strengthening occurs within the first few seconds of the application of the restraining force and is highly localized.Cells ‘sense’ the stiffness of their substrates and exert smaller tractions on more compliant substrates (Lo

*et al*. 2000; Discher*et al*. 2005).

## 3. Constitutive model for a single stress fibre

### (a) The model

The preceding biochemistry suggests that the mechanical response of the stress fibres comprises three coupled phenomena:

An activation signal that triggers the formation of actin and bipolar myosin filaments.

A fibre-formation rate dependent on the activation signal, coupled with a dissociation rate affected by the tension.

A contraction rate (contractility) for the stress fibre that depends on the tension through the cross-bridge dynamics.

Phenomenological relations are proposed for each of these features. The expressions chosen are the simplest possible without violating the biochemistry. In future developments, when warranted, different mathematical dependencies can be considered for all three. These will not alter the general features predicted by the model: only the absolute magnitudes and the time-scales will change.

#### (i) Rapid transmission of the extracellular signal triggering the polymerization of the actin filaments and the phosphorylation of the myosin

The precise details of the signalling pathways are ignored. Rather, the level, which may be thought of as the concentration of Ca^{2+}, designated *C*(0≤*C*≤1), is assumed given by(3.1)where *θ* is the decay constant of the signal and *t*_{i} is the time measured from the instant of the *most recent* signal. Over the time-scales of the contractile activity of the cell (of the order of hours), we assume that there are no spatial gradients in *C* and that the level is dominated by the most recent signal. Readers are referred to Rudy (2000) for more sophisticated models of the Ca^{2+} dynamics.

#### (ii) Polymerization/depolymerization of the actin filaments and phosphorylation/dephosphorylation of the myosin

The transduction of the signal results in: (i) the polymerization of the actin filaments and the bundling of these filaments by α-actinin and (ii) the phosphorylation of myosin II, which promotes the assembly of the myosin into bipolar filaments. The interaction between the myosin II heads and the actin filaments forms contractile bundles.

We characterize the activation level of the stress fibre bundles by a non-dimensional parameter, *η* (0≤*η*≤1), defined as the ratio of the concentration of the polymerized actin and phosphorylated myosin in the bundle to the maximum concentrations permitted by the biochemistry (Deshpande *et al*. 2006). The formation and the dissociation of the stress fibres, as parameterized through *η*, are represented by a first-order kinetic scheme,(3.2)

The over-dot denotes differentiation with respect to time, *t*, measured from the instant of the application of the first signal. The term *T* is the tension in the stress fibre and *T*_{0}(*η*) the corresponding isometric tension for a given *η*. The dimensionless constants and govern the rate of formation and dissociation of the stress fibres, respectively. The key features of equation (3.2) are as follows: (i) the rate of the formation of the stress fibre (first square bracket) decreases with increasing fibre activation *η* and is proportional to the strength of the decaying signal and (ii) the rate of dissociation (second square bracket) is proportional to the concentration of the polymerized actin and phosphorylated myosin II. Moreover, to allow the contractile bundles to be held together by the tension, we propose a dissociation term dependent on the tension; namely, the dissociation rate is zero when the fibres are held at their isometric tension, but increases linearly at lower tension.

The isometric tension *T*_{0} is taken to be directly proportional to the level of activation of the stress fibre, *T*_{0}=*ηT*_{max}. Here, *T*_{max} is the isometric tension when the concentration of the polymerized actin and phosphorylated myosin is the maximum permitted by the biochemistry. We note that, in muscle cells, the isometric tension depends on the overlap between the thick and the thin filaments giving rise to a dependence on the muscle length. By contrast, in non-muscle cells, cytoskeletal rearrangements, along with the formation of focal adhesions, govern the dynamic formation and dissociation of the stress fibres; namely, given the absence of a static cytoskeletal structure, there is no physical motivation to include length dependence.

#### (iii) Tension versus velocity relation and the cross-bridge dynamics

The tension in the stress fibres is generated by the cross-bridge cycling between the actin and the myosin filaments. This force generation mechanism is similar (but not identical) to that in muscle cells. Consequently, we assume that the influence of tension on the extension/shortening rate of the fibres is adequately described by a Hill-like relation (Hill 1938)2. The Hill measurements (sketched in figure 2) motivate a simple tension versus velocity relationship,(3.3)Here, *v* is the rate of change of the length of stress fibre (positive for lengthening and negative for shortening). The non-dimensional constant, , is the fractional reduction in tension when the shortening rate increases by the reference value, *v*_{0}. Under shortening conditions, this relation captures the main features of the Hill relation (figure 2). Alternatives to equation (3.3) are discussed in Appendix B in the electronic supplementary material.

A passive elastic contribution could also be included in the one-dimensional stress fibre relation by assuming that the total force is the sum of the active force as given by equations (3.1)–(3.3) and a force due to elastic stretching of elements such as the intermediate filaments that act in parallel to the stress fibre. Such passive elastic contributions will be considered in the two- and three-dimensional models presented subsequently, but are neglected in the one-dimensional model in order to better illustrate the essential features of the stress fibre constitutive relation.

### (b) Illustrative examples

To give physical insight into these coupled equations, some examples are presented for a single stress fibre bundle (with full realization that these examples do not represent realistic situations). For simplicity, we assume that the fibre is completely inactive at time *t*=0. Thus, the response reduces to solving the ordinary differential equations (3.1)–(3.3) for given boundary conditions with the initial condition *η*=0 at *t*=0.

*Effect of support stiffness*. Consider a stress fibre held between two rigid foundations through a support spring in series (figure 3), characterized by a spring constant *k*_{s}. This spring constant is defined via the relation *T*=*k*_{s}Δ*x*, where *T* is the tension generated by the spring for an extension Δ*x*. (Force equilibrium dictates equality of the tension in the stress fibre and the spring.) Activation initiates shortening of the stress fibre resisted by the support springs. The resulting coupled governing equations are equation (3.2) and(3.4)with initial conditions *η*=*T*=0 at *t*=0. The problem is fully specified by four independent non-dimensional variables: , , the non-dimensional support stiffness,(3.5)and the non-dimensional time, . The resulting solution-dependent variables are the non-dimensional tension in the stress fibre (and spring), , and the non-dimensional extension of the support spring, . Unless otherwise specified, the reaction rate constants are taken as and , referred to as the reference case. The non-dimensional spring constant is varied in the range .

The effect of the spring stiffness on the evolution of tension with time is plotted in figure 4*a*. For stiff springs (), *the steady-state tension always reaches T*_{max}. Conversely, compliant springs () are unable to sustain significant tension until the stress fibre has undergone substantial contraction. Whereupon, *η* increases more slowly and the stress fibre does not achieve full activation before the signal *C* has decayed. Consequently, *the steady-state levels of* *η* *and* *are both less than unity*. The variation of the steady-state tension with (figure 4*b*) reveals the interplay between and the two reaction rates, and . Increasing the forward rate delays the drop-off in to lower . Increasing the dissociation rate of the fibres causes to fall below its maximum at higher .

*A stress fibre on an array of N posts*. This problem is motivated by experimental investigations of the responses of cells on two-dimensional arrays of micro-needles (Tan *et al*. 2003). These investigations revealed that the average force per post increases with increasing number of posts *N*. Here, we consider the analogous one-dimensional problem (figure 5) of a stress fibre spanning *N* equispaced posts, each with stiffness *k*_{s} (force *F*_{i} on the *i*th post is related to its displacement via *F*_{i}=*k*_{s}Δ*x*_{i}, where Δ*x*_{i} is the displacement of the *i*th post). The governing equations are derived here for a stress fibre on an odd number, *N*=2*M*+1, of posts (an analogous set of equations follow for an even number of posts). The array is symmetrical about the central post, labelled ‘post 0’ and it suffices to consider only one-half of the array (the posts are labelled as indicated in figure 5). With the coordinate *x* increasing from left to right (figure 5), the equilibrium and the compatibility equations for the right half of the array are given, respectively, as(3.6a)(3.6b)where *T*_{i} and *v*_{i} are the tension and the extension rate of the *i*th fibre segment between the (*i*−1)th and *i*th posts. These equilibrium and compatibility conditions are combined with the constitutive equations(3.7a)and(3.7b)and solved using a fifth-order Runge–Kutta method with initial conditions *η*_{i}=*T*_{i}*=*Δ*x*_{i}*=0* at *t*=0 (*η*_{i} and *v*_{i} are the activation level and extension rate, respectively, of the *i*th fibre segment). The forward and backward reaction rate constants are fixed at and we investigate the effects of *N* and the post stiffness *k*_{s}. With defined as the stead-state force on the *i*th post, the non-dimensional average steady-state force per post is(3.8)

Plots for three choices of the non-dimensional post stiffness (figure 6) reveal that first increases with *N*, reaches a maximum and then decreases. The maximum is attained at larger *N* with decreasing . The origin of this scaling emerges from the following considerations. For a stress fibre spanning three posts, the contraction is resisted only by the two outer posts (by symmetry the central post does not deflect). For larger *N*, the contraction of the fibre segments adjacent to the centre is resisted by an increasing number of posts (equation (3.6*a*)), enhancing the resistance to contraction and, thereby, elevating (figure 4). With still more posts, all interior fibre segments achieve their maximum tension *T*_{max}, eliminating the force gradients towards the mid-section. Consequently, the forces on the central posts approach zero, causing a reduction in the net force with increasing *N*.

To further clarify, the temporal dependence of the deflections, stresses and forces are examined. The post deflections for the 11-post case (figure 7) reveal that, in steady-state, the outermost post undergoes the largest deflections, while those at the inner posts become increasingly smaller, corresponding to the edge shear-lag phenomenon identified by Mohrdieck *et al*. (2005) and discussed above. The associated effects are revealed by plotting the lengthwise distribution of the steady-state tension in the stress fibre (figure 8*a*), and the corresponding non-dimensional steady-state post forces, , where is the steady-state post deflection (figure 8*b*). The gradients in the tension in the fibre first increase with increasing *N* (up to *N*=5) but then start to reduce towards the interior. This tension distribution in the fibre results in the forces on the posts near the centre decreasing with increasing *N* and nearly reducing to zero when *N*>5.

## 4. A constitutive model for a two-dimensional cytoskeletal network

We generalize the model to a cell comprising a two-dimensional network of stress fibres by invoking the following key assumptions.

There is sufficient actin and myosin in the cell that the activation of the stress fibres in each direction is not limited by their availability. Rather, it is limited by the tension and polymerization, expressed through equation (3.2).

A representative volume element (RVE)3 on a length-scale much smaller than any leading dimension of the cell can be defined (figure 9). The RVE has pillbox shape with radius

*R*, depth*h*and volume*V*≡*πR*^{2}*h*.Stress fibres can form in any direction

*ϕ*with equal probability.

By adopting Cartesian tensor notation (Einstein summation convention over repeated indices), the average Cauchy stress *σ*_{ij} over the RVE resulting from the fibres is defined from the stress *τ*_{ij} at any point as(4.1)with respect to the fixed orthogonal basis *x*_{i} (figure 9). This RVE contains a large number of stress fibres at different levels of polymerization, uniformly distributed over −*π*/2≤*ϕ*≤*π*/2, where *ϕ* is the fibre angle measured with respect to an orthogonal set of base vectors *e*_{i} that rotate (but do not deform) with the material. Since *e*_{i} are coincident with at time *t*=0, *ϕ* gives the stress fibre orientation in the original configuration. Unless otherwise stated, all tensor and vector components are measured with respect to the fixed basis .

The virtual work statement for the RVE is(4.2)where is the volume average of the strain rate (the symmetric part of the spatial velocity gradient), is an arbitrary variation in and *δv* is the corresponding variation in the fibre extension rate, *v*. Employing compatibility, the axial strain rate in the fibre, , is related to by(4.3)Here, *ϕ*^{*} is the angle of the stress fibre measured with respect to and is related to *ϕ* by the rigid body rotation of the material point. In a two-dimensional setting (with anticlockwise rotations taken to be positive),(4.4)where are the displacement rates. Upon combining equations (4.1) and (4.3), the virtual work statement (4.2) reduces to(4.5)Recall that are arbitrary variations in . Thus, by successively allowing one component of to be non-zero while the others are set to zero (e.g. to obtain *σ*_{11}), the average stress is given as(4.6a)where the ‘smeared-out’ stress *σ* is related to the fibre tension *T* via(4.6b)We now prescribe the governing equations for *σ*, analogous to the evolution equations for *T* for a single stress fibre. Recall that the activation of the stress fibres in each direction is regarded as decoupled and not limited by the availability of actin or myosin. Thus, the evolution of the activation level at an angle *ϕ* becomes(4.7a)where the isometric stress is *σ*_{0}(*ϕ*)=*η*(*ϕ*)*σ*_{max}. Here, *σ*_{0} and *σ*_{max} are related to the fibre tensions via(4.7b)respectively. With the fibre extension and strain rate related by , the Hill-like equation (3.3) is generalized to the intensive quantities of the fibre stress *σ* and fibre strain rate as(4.8)The active properties of the stress fibres, as parameterized through *σ*_{max}, , and , remain unaffected by the material deformation as the cytoskeletal rearrangements essentially ‘reset’ the stress fibre properties with continued deformation. Consequently, a ‘natural’ state for the stress fibres does not exist.

The contractile response in two dimensions is expected to include a contribution from the passive elasticity provided mainly by the intermediate filaments of the cytoskeleton attached to the nuclear and plasma membranes. We assume additive decomposition of the active and passive stresses since the fibres act in parallel with the intermediate filaments and the cell membrane. Here, for simplicity, we employ a finite strain elasticity that readily reduces to the isotropic linear elastic Hooke's law for infinitesimal deformation. The passive elastic response is described by a strain energy density function *W*_{0} (per unit undeformed volume),(4.9a)where *E* and *ν* are Young's modulus and Poisson's ratio, respectively, and is the Green–Lagrange strain(4.9b)Here, *δ*_{ij} is the Kronecker delta and the deformation gradient is given as *F*_{ij}=*δ*_{ij}+(∂*u*_{i}/∂*x*_{j}). The work conjugate to the Green–Lagrange strain is the second Piola–Kirchhoff stress , where *Σ*_{ij} is the Cauchy stress and *J*=det (*F*_{ij}). Thus, the passive elastic stress contribution is specified via(4.10)Note that, in the numerical examples presented here and in Deshpande *et al*. (2006), *the strains in the cell are relatively small and a linear elastic relation for the passive elasticity suffices*. When warranted, equation (4.10) can readily be replaced by a nonlinear (hyperelastic) law. Employing additive decomposition, the total Cauchy stress *Σ*_{ij}, combining active contributions from equation (4.6*a*) and passive elasticity from equation (4.10), is obtained by transforming , such that(4.11)The three-dimensional analogue of the active features of the constitutive law is presented in Appendix C (see electronic supplementary material), for which equation (4.11) can then be used as given.

The foregoing model has been implemented as a user-defined material model (UMAT) in the commercial finite element package ABAQUS4. Details of the numerical implementation are given in Appendix D in the electronic supplementary material. All calculations are performed in plane stress (*σ*_{33}=*Σ*_{33}=0) finite deformation setting, i.e. the effects of geometry changes on the momentum balance and rigid body rotations are taken into account. Four-noded elements (CPS4 in ABAQUS notation) are employed.

### (a) Reference material properties

No attempt is made to justify the choice of the parameters employed using either theoretical arguments or precise experimental measurements. Rather, these constants have been chosen to give results similar to those in Tan *et al*. (2003). The decay constant of the signal was taken to be *θ*=720 s, while the passive Young's modulus and Poisson's ratio chosen to be *E*=0.077 nN μm^{−2} and *ν*=0.3, respectively. It is worth emphasizing here that the permeability of the cell walls implies that cells typically do not exhibit an incompressible response. Hence, unlike soft tissues, the choice *ν*≈0.5 is not appropriate for single cells. The non-dimensional reaction rate constants are and , while the non-dimensional fibre rate sensitivity is . The maximum tension exerted by the stress fibres is *σ*_{max}=3.9 nN μm^{−2} and the reference strain rate in the cross-bridge dynamics law is . The non-dimensional material properties of the cell are . Unless otherwise specified, in the results presented subsequently, the corresponding non-dimensional cell properties are , *ν*=0.3, , , and . The cell is initially stress and stress-fibre free; namely, the initial conditions are *η*(*ϕ*)=0 at time *t*=0 over the entire cell.

### (b) Development of structural anisotropy in a cell under isometric conditions

Here, we demonstrate the capabilities of the constitutive model for predicting the development of structural anisotropy. Rectangular cells, dimensions *l*×*h* (figure 10), are subjected to a single activation signal at time *t*=0 with the following boundary conditions: (i) *Uniaxial isometric tension*: *u*_{1}=0 on the left and the right edges of length *h*, with tractions *T*_{1}=*T*_{2}=0 on the top and the bottom surfaces of length *l* and *T*_{2}=0 on the left and the right edges (figure 10*a*). (ii) *Biaxial isometric tension*: *u*_{1}=0 on the left and the right edges of length *h* and *u*_{2}=0 on the top and the bottom surfaces of length *l*. The tractions, *T*_{2}=0 and *T*_{1}=0, on the vertical and horizontal edges, respectively (figure 10*b*). Under both loading conditions, the fields () are spatially uniform with no rotation of the material. All results are thus independent of the dimensions of the cell.

For *uniaxial isometric loading*, the time evolution of the normalized stress *Σ*_{11}/*σ*_{max} (figure 11*a*) and the corresponding evolution of the logarithmic transverse strain (figure 11*b*) indicate that, after an initial peak stress at *t*/*θ*≈1, a steady state is achieved at *t*/*θ*≈30, coincident with cessation of the transverse straining of the cells. The corresponding distributions of stress fibre activation levels *η*(*ϕ*) are shown in figure 12*a* at four selected *t*/*θ* (the steady-state distribution of *η* is marked as *t*/*θ*=∞). Recall that the material does not undergo rigid body rotation, so that the stress fibre orientations *ϕ* in the undeformed configuration are the same as their respective orientations *ϕ** in the deformed configuration. At steady state, the stress fibres tend to align with the loading direction, resulting in a highly anisotropic cytoskeletal structure (figure 12*b*). The results can be understood as follows. Immediately after the signal is applied, stress fibres are activated in all directions. Fibres at *ϕ*≈0° are constrained against shortening by the uniaxial isometric boundary condition (albeit to different extents). They remain activated, as evident in figure 12*a* at *t*/*θ*=0.5. However, fibres at *ϕ*≈90° are unconstrained and shorten at strain rate , with negligible *σ*_{22} stress (as required by the traction-free boundary conditions on the top and the bottom surfaces of the cell). Thus, the cell starts to contract in the *x*_{2}-direction. But since a small stress *σ*_{22} is generated in *x*_{2}, these fibres begin to dissociate. The outcome at steady state is as follows. The fibres at *ϕ*=0° are fully activated (*η*=1.0), while those at |*ϕ*|<20° remain active, but with 0<*η*<1.0, and all others completely dissociate (*η*=0). The tensile *σ*_{22} stresses generated by the fibres at *ϕ*≠0° are balanced by the elastic stresses due to compression in the *x*_{2}-direction. Increasing elevates the isomeric tension and decreases the transverse straining (figures 11 and 12). Moreover, with increasing steady state is achieved sooner, without overshoot in the generated tension. This trend arises because the fibres orthogonal to the tensile axis contract and generate an axial compressive stress (owing to the Poisson effect), thereby reducing the axial tension. Increasing Young's modulus decreases this contraction in the *x*_{2}-direction. The effect of on the structural anisotropy is illustrated in figure 12, revealing that increasing causes fibres over a wider range of *ϕ* to be recruited into the tension generating apparatus, thereupon reducing the structural anisotropy.

*For the cell under biaxial loading*, since no straining occurs, the passive Young's modulus does not affect the stresses generated. The cell generates the highest tension at *t*/*θ*≈1 (figure 11*a*) with a steady-state stress, *Σ*_{11}/*σ*_{max}=*Σ*_{22}/*σ*_{max}=0.5.

## 5. Responses of a square cell on four supports

Experiments to probe the forces generated by a cell on a bed of micro-needles (Tan *et al*. 2003) have motivated the two-dimensional plane stress problem illustrated in figure 13. Using this representation, the focal adhesions need not be considered explicitly. A similar problem considered by Deshpande *et al*. (2006; figures 4–7) demonstrated that the foregoing constitutive model captures key experimental observations: (i) the decrease of the forces generated by the cell with decreasing support stiffness, (ii) the influence of cell shape and boundary conditions on the development of structural anisotropy, and (iii) the high concentration of the stress fibres at the focal adhesions. Here, we demonstrate two additional features: (i) the effect of multiple activation signals and (ii) the response to a localized force.

### (a) Geometry and boundary conditions

A square cell with reference properties, side *L*=50 μm (thickness *b*=1 μm), is supported over a length *L*_{s}=5 μm at the four edges by identical elastic foundations of stiffness *k*_{E}=0.77 nN μm^{−3} (figure 13). The foundation rotates with the cell edges and thus can only exert traction normal to the supported cell edges, i.e. in this finite deformation setting, the traction rate on a cell surface with unit outward normal (in the current configuration) is given by(5.1)where is the displacement rate along and is the magnitude of the traction vector. Initially, (*t*=0), *η*(*ϕ*)=0 over the entire cell, with the first activation signal applied at *t*=0. A uniform FE mesh with an element size of 0.25 μm was employed in all calculations. The additional non-dimensional variables in the problem are and . In all the examples presented here, .

Results are presented for the time evolution of the average displacement over the support area, defined as(5.2)where *u*_{i} are the displacement components of the cell along the supported edges. The corresponding work-conjugate non-dimensional support force is . In order to visualize the evolution of the stress fibres, we define two additional quantities.

The average stress fibre activation over all orientations,(5.3)

The maximum principal value of the active stresses,

*σ*_{ij}, and the associated principal direction, measured as the orientation*ϕ*_{p}with respect to the*x*_{1}-axis (figure 13). The orientation*ϕ*_{p}may be regarded as the ‘resultant’ stress fibre direction.

### (b) Effect of multiple activations

The motivation is the suggestion that a single long activation of osteoblasts (cells that play a critical role in the process of new bone formation) is less successful in reorganizing the cytoskeleton into stress fibres than multiple short activations separated by rest periods (Robling *et al*. 2001). This suggestion has been made based on the observation that bone formation is enhanced when activations are interspersed with rest periods. Here, we contrast the multiple and single activation responses of the cell.

In order to contrast the response of the cell subjected to a single sustained and multiple activations, we consider the following three activation prescriptions for the cell of figure 13*a*.

A single activation signal with a decay constant of

*θ*=1440 s applied to the cell at time (i.e. with all other cell properties fixed at their reference values).Two activation signals (both with a decay constant of

*θ*=720 s) applied to the cell at times . All the cell properties are thus unchanged from their reference values in this case.Four activation signals (each with a decay constant of

*θ*=360 s) applied to the cell at times , 15, 30 and 45. Thus, with all other cell properties fixed at their reference values.

These three prescriptions have been chosen such that the total activation time *Σθ*=1440 s in all the cases and in the case of multiple activations, the next activation is only applied after the response of the cell has attained a steady state as a result of the previous activation.

The time evolutions of the average support displacement, , and forces, (figure 14), indicate that while the final steady state attained by the cell is similar for two and four activation signals, these steady-state levels are higher than those attained by the cell when it is subjected to a single prolonged activation. Moreover, when the cell is subjected to a single prolonged activation, the force exerted by the cell at the supports initially overshoots, before decreasing towards the steady state. The distributions of for the single prolonged activation (*θ*=1440 s) and the two activation signal (*θ*=720 s) cases are plotted in figure 15 (at *t*/*θ*=15 and 60) with the orientations of the principal stress included as line segments whose length is scaled by the normalized stress, . The long activation signal results in a higher concentration of stress fibres at , while at steady state (), higher concentrations develop in the cell subjected to two activation signals. This indicates that the nonlinearity of the cell response results in the cell rearranging the cytoskeleton more effectively when subjected to multiple activations rather than one prolonged activation.

### (c) Response of the cell to the application of a localized force

Experiments by Kolega (1986) have demonstrated that tension applied to a localized site on the cell surface results in a bundle of stress fibres being induced immediately adjacent to this site. Motivated by this observation, we consider the problem sketched in figure 13*b* where a localized tension is imposed over a central patch on the right-hand surface of the cell, in accordance with the following conditions: (i) a single activation signal is applied at *t*=0 and (ii) a displacement rate is applied in the *x*_{1}-direction over a centrally situated patch, length *L*_{f}=5 μm, on the right surface of the cell (figure 13*b*) over duration, (where ). Subsequently, the patch is fixed, with displacement rate prescribed in the *x*_{1}-direction, over *L*_{f}.

The work-conjugate force associated with the applied displacement rate is designated *F*_{f}. The following additional non-dimensional groups are present here: (i) the patch length , (ii) the loading rate , and (iii) the force that is the work conjugate to the applied displacement rate . In the calculations, with displacement rate applied over the duration and for .

The normalized force required to continue displacing the patch at a rate increases with time (figure 16). When the displacement rate is set to zero (at ), the force begins to decrease before reaching steady state, , at time . The distributions of the average stress fibre activation levels at (figure 17*a*,*b*, respectively) with the orientations of reveal that (i) a high concentration of stress fibres is induced immediately adjacent to the site of applied tension, (ii) stress fibres form normal to the surface at this location, and (iii) at steady state (), a high concentration of stress fibres continues to persist near the supports and at the site of the applied tension.

## 6. Discussion

Contractility in non-muscle cells has previously been modelled by prescribing a thermal strain to either a cell otherwise regarded as an isotropic elastic continuum (Nelson *et al*. 2005) or a discrete set of elastic filaments representing the stress fibres (Mohrdieck *et al*. 2005). Such models neglect the biochemistry of the active apparatus of the cell that generates, supports and responds to mechanical forces. The model presented here accounts for the nonlinear coupling between signalling, the kinetics of tension-dependent stress-fibre formation/dissolution and stress-dependent contractility. The distinguishing features of this model are as follows.

It accounts for micro-structural evolution of the cytoskeleton. In particular, the dynamic reorganization of the cytoskeleton into stress fibres has been motivated through biochemical considerations coupled with the mechanics.

The active apparatus of a cell cannot support compressive forces or stresses and thus Ingber (1997) suggested ‘tensegrity’ as the architectural basis for the cytoskeletal arrangements. The current model automatically ensures that only tensile active stresses are generated without assuming

*a priori*an arrangement of stress fibres that ensures a tensegrity structure. Rather, the tensile stress fibre network is allowed to evolve in a manner dictated by the boundary conditions applied to the cell.

These critical features enable the model to predict a wide range of experimentally observed phenomena. However, a note of caution is appropriate. The model is highly nonlinear and the results depend on the choice of the parameters. For example, the results in figure 14 indicate that the steady-state forces due to two or four activations are reasonably similar, while a single prolonged activation results in significantly smaller support forces. This outcome will depend strongly on the couplings between the signalling, fibre formation kinetics, dissociation and contractility. Large-scale parametric studies will be used to fully understand these nonlinear couplings before embarking on generalized statements.

## 7. Concluding remarks

A biochemically inspired model has been presented for the dynamic rearrangement of the cytoskeleton that incorporates cell contractility. The constitutive equations are presented for the formation and response of a single stress fibre and a model developed for a two- and three-dimensional cytoskeletal network. The contribution of passive elasticity is also included in the two- and three-dimensional models. All models entail the highly nonlinear interaction between signalling, the kinetics of tension-dependent stress-fibre formation/dissolution and stress-dependent contractility.

One-dimensional numerical examples of a single stress fibre constrained between rigid supports with a spring in series, as well as a stress fibre on an array of posts, are presented in order to illustrate the salient features of the model including: (i) the force exerted by the stress fibre increases with increasing support stiffness and (ii) the average force exerted by the stress fibre on the array of post first increases with increasing number of posts and thereafter decreases. Consistent with experimental observations, the maximum deflection of the posts occurs towards the edges of the array.

Numerical examples have been presented for the contractility of a cell with a two-dimensional cytoskeletal network. The uniaxial and biaxial isometric tension responses have been investigated. The model captures the development of structural anisotropy under uniaxial isometric loading conditions, while the cell remains isotropic under biaxial loading. It predicts a strong coupling between the development of anisotropy and the passive elasticity of the cell: increasing the passive Young's modulus diffuses the formation of stress fibres and reduces the anisotropy.

The response of a square cell on four corner supports has been considered. In line with the experimental findings, the model predicts that (i) a single prolonged activation signal is less effective at developing stress fibres than multiple shorter signals and (ii) a high concentration of stress fibres is formed near the supports and the sites of localized applied tension.

The results taken together with those in Deshpande *et al*. (2006) demonstrate that this model captures a wide range of key experimental findings on the contractility of non-muscle cells. The correspondences between the simulations presented in Deshpande *et al*. (2006) and the *in vitro* observations include: (i) the decrease of forces generated by the cell with increasing substrate compliance, (ii) the influence of cell shape and boundary conditions on the development of structural anisotropy, and (iii) the high concentration of stress fibres at focal adhesions. However, direct contact with experiments such as those in Parker *et al*. (2002) and Tan *et al*. (2003) will require (i) calibration of the constants in the model against experimental data and (ii) the inclusion of additional features in the model such as the effects of cell spreading and the development of focal adhesions.

The model can be used to address one of the key challenges in cell biomechanics; namely, how to measure the mechanical characteristics of living cells that react to the measurement tools. Since the model captures the reorganization of the cytoskeletal elements in response to mechanical perturbations, it can be employed as a framework to design and interpret appropriate experiments.

## Acknowledgments

V.S.D. acknowledges support from the Leverhulme Trust, UK. R.M.M. and A.G.E. thank the Army Research Office for their support through a subcontact with the University of Chicago provided by a Multidisciplinary University Research Initiative Program on ‘Bio-Mechanical Interfaces for Cell-Based Microsystems’, Prime Award No.: W911NF-04-1-071. The authors are grateful to Professors C. Chen, A. Dinner and M. Mrksich for their helpful discussions and comments.

## Footnotes

Electronic supplementary material is available at http://dx.doi.org/10.1098/rspa.2006.1793 or via http://www.journals.royalsoc.ac.uk.

↵In contrast to non-muscle cells, the cytoskeleton in muscle cells is relatively static and comprises a regular array of filament bundles called sarcomeres.

↵The Hill (1938) equation adequately models the cross-bridge dynamics in the low-frequency domain and readers are referred to Rice

*et al*. (1999) for a discussion on more sophisticated models that capture the cross-bridge dynamics over a large frequency range. Since cell contractility is a relatively slow process, it suffices to use the simple Hill-type description to model the tension versus velocity relationship of the stress fibres.↵Hill (1963) defined an RVE as a domain that is (i) structurally typical of the entire mixture (network of stress fibres) and (ii) sufficiently large compared with the micro-scale (say, the stress fibre radius).

↵ABAQUS (2004) User's Manual, v. 6.5. ABAQUS, Inc.

- Received May 18, 2006.
- Accepted November 14, 2006.

- © 2007 The Royal Society