A model for the contractility of the cytoskeleton including the effects of stress-fibre formation and dissociation

Vikram S Deshpande, Robert M McMeeking, Anthony G Evans

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.

Keywords:

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.

Figure 1

Measurement of contractile forces in a fibroblast cell on a bed of micro-needles. The actin fibres are stained in green. The arrows show the deflection of the posts with the lengths of the arrows proportional to the force exerted by the cell on the posts. There seems little correlation between the orientations of the visible stress fibres and the directions of the force vectors (courtesy of Prof. C. Chen, University of Pennsylvania).

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.

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

  2. 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).

  3. 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 PIP2 and stimulates the release of Ca2+ from the endoplasmic reticulum. The influx of Ca2+ 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 Ca2+ 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

  1. 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).

  2. 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.

  3. 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:

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

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

  3. 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 Ca2+, designated C(0≤C≤1), is assumed given byEmbedded Image(3.1)where θ is the decay constant of the signal and ti 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 Ca2+ 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,Embedded Image(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 T0(η) the corresponding isometric tension for a given η. The dimensionless constants Embedded Image and Embedded Image 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 T0 is taken to be directly proportional to the level of activation of the stress fibre, T0=ηTmax. Here, Tmax 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,Embedded Image(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, Embedded Image, is the fractional reduction in tension when the shortening rate increases by the reference value, v0. 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.

Figure 2

Hill (1938) tension, velocity relation for skeletal muscles. The simple approximation (3.3) to the Hill relation used in this study is also included in the figure for comparison purposes.

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 ks. This spring constant is defined via the relation T=ksΔ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) andEmbedded Image(3.4)with initial conditions η=T=0 at t=0. The problem is fully specified by four independent non-dimensional variables: Embedded Image, Embedded Image, the non-dimensional support stiffness,Embedded Image(3.5)and the non-dimensional time, Embedded Image. The resulting solution-dependent variables are the non-dimensional tension in the stress fibre (and spring), Embedded Image, and the non-dimensional extension of the support spring, Embedded Image. Unless otherwise specified, the reaction rate constants are taken as Embedded Image and Embedded Image, referred to as the reference case. The non-dimensional spring constant is varied in the range Embedded Image.

Figure 3

An isolated stress fibre fixed between two rigid supports via a spring in series.

The effect of the spring stiffness on the evolution of tension Embedded Image with time Embedded Image is plotted in figure 4a. For stiff springs (Embedded Image), the steady-state tension always reaches Tmax. Conversely, compliant springs (Embedded Image) 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 Embedded Image are both less than unity. The variation of the steady-state tension Embedded Image with Embedded Image (figure 4b) reveals the interplay between Embedded Image and the two reaction rates, Embedded Image and Embedded Image. Increasing the forward rate delays the drop-off in Embedded Image to lower Embedded Image. Increasing the dissociation rate of the fibres causes Embedded Image to fall below its maximum at higher Embedded Image.

Figure 4

Effect of the spring stiffness on the response of the stress fibre. (a) The time evolution of the normalized tension Embedded Image with non-dimensional time Embedded Image for four choices of the support stiffness Embedded Image, using the reaction constants: Embedded Image and Embedded Image. (b) The normalized steady-state tension, Embedded Image, as a function of Embedded Image for three combinations of the forward and backward rate constants Embedded Image and Embedded Image.

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 ks (force Fi on the ith post is related to its displacement via Fi=ksΔxi, where Δxi is the displacement of the ith post). The governing equations are derived here for a stress fibre on an odd number, N=2M+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, asEmbedded Image(3.6a)Embedded Image(3.6b)where Ti and vi are the tension and the extension rate of the ith fibre segment between the (i−1)th and ith posts. These equilibrium and compatibility conditions are combined with the constitutive equationsEmbedded Image(3.7a)andEmbedded Image(3.7b)and solved using a fifth-order Runge–Kutta method with initial conditions ηi=Ti=Δxi=0 at t=0 (ηi and vi are the activation level and extension rate, respectively, of the ith fibre segment). The forward and backward reaction rate constants are fixed at Embedded Image and we investigate the effects of N and the post stiffness ks. With Embedded Image defined as the stead-state force on the ith post, the non-dimensional average steady-state force per post isEmbedded Image(3.8)

Figure 5

One-dimensional array of identical posts spaced a distance l apart with a stress fibre spanning these posts. A seven-post array is shown here with the central post designated post 0. The ith fibre segment resides between the ith and (i−1)th posts.

Plots for three choices of the non-dimensional post stiffness (figure 6) reveal that Embedded Image first increases with N, reaches a maximum and then decreases. The maximum Embedded Image is attained at larger N with decreasing Embedded Image. 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.6a)), enhancing the resistance to contraction and, thereby, elevating Embedded Image (figure 4). With still more posts, all interior fibre segments achieve their maximum tension Tmax, eliminating the force gradients towards the mid-section. Consequently, the forces on the central posts approach zero, causing a reduction in the net force Embedded Image with increasing N.

Figure 6

Non-dimensional steady-state average force per post Embedded Image as a function of the number of posts N for three choices of the post spring constant Embedded Image. The forward and backward rate constants of the stress fibre are Embedded Image.

To further clarify, the temporal dependence of the deflections, stresses and forces are examined. The post deflections Embedded Image 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 Embedded Image in the stress fibre (figure 8a), and the corresponding non-dimensional steady-state post forces, Embedded Image, where Embedded Image is the steady-state post deflection (figure 8b). 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.

Figure 7

Normalized post deflection versus time history for a stress fibre on an array of 11 posts. The forward and backward rate constants of the stress fibre are Embedded Image and the post spring stiffness is Embedded Image.

Figure 8

(a) Distribution of the normalized steady-state tension Embedded Image along the stress fibre for 3-, 5-, 7- and 11-post arrays and (b) steady-state force Embedded Image on each post for the stress fibre on 3-, 5- and 11-post arrays. In both plots, the forward and backward rate constants are Embedded Image and the post spring stiffness is Embedded Image.

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.

  1. 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).

  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πR2h.

    Figure 9

    Macro- and micro-scales in a cell with a two-dimensional network of stress fibres. The representative volume element (RVE) in the shape of a pillbox is also shown.

  3. 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 asEmbedded Image(4.1)with respect to the fixed orthogonal basis xi (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 ei that rotate (but do not deform) with the material. Since ei are coincident with Embedded Image 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 Embedded Image.

The virtual work statement for the RVE isEmbedded Image(4.2)where Embedded Image is the volume average of the strain rate (the symmetric part of the spatial velocity gradient), Embedded Image is an arbitrary variation in Embedded Image and δv is the corresponding variation in the fibre extension rate, v. Employing compatibility, the axial strain rate in the fibre, Embedded Image, is related to Embedded Image byEmbedded Image(4.3)Here, ϕ* is the angle of the stress fibre measured with respect to Embedded Image 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),Embedded Image(4.4)where Embedded Image are the displacement rates. Upon combining equations (4.1) and (4.3), the virtual work statement (4.2) reduces toEmbedded Image(4.5)Recall that Embedded Image are arbitrary variations in Embedded Image. Thus, by successively allowing one component of Embedded Image to be non-zero while the others are set to zero (e.g. Embedded Image to obtain σ11), the average stress is given asEmbedded Image(4.6a)where the ‘smeared-out’ stress σ is related to the fibre tension T viaEmbedded Image(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 ϕ becomesEmbedded Image(4.7a)where the isometric stress is σ0(ϕ)=η(ϕ)σmax. Here, σ0 and σmax are related to the fibre tensions viaEmbedded Image(4.7b)respectively. With the fibre extension and strain rate related by Embedded Image, the Hill-like equation (3.3) is generalized to the intensive quantities of the fibre stress σ and fibre strain rate Embedded Image asEmbedded Image(4.8)The active properties of the stress fibres, as parameterized through σmax, Embedded Image, Embedded Image and Embedded Image, 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 W0 (per unit undeformed volume),Embedded Image(4.9a)where E and ν are Young's modulus and Poisson's ratio, respectively, and Embedded Image is the Green–Lagrange strainEmbedded Image(4.9b)Here, δij is the Kronecker delta and the deformation gradient is given as Fij=δij+(∂ui/∂xj). The work conjugate to the Green–Lagrange strain is the second Piola–Kirchhoff stress Embedded Image, where Σij is the Cauchy stress and J=det (Fij). Thus, the passive elastic stress contribution is specified viaEmbedded Image(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.6a) and passive elasticity from equation (4.10), is obtained by transforming Embedded Image, such thatEmbedded Image(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 Embedded Image and Embedded Image, while the non-dimensional fibre rate sensitivity is Embedded Image. 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 Embedded Image. The non-dimensional material properties of the cell are Embedded Image. Unless otherwise specified, in the results presented subsequently, the corresponding non-dimensional cell properties are Embedded Image, ν=0.3, Embedded Image, Embedded Image, Embedded Image and Embedded Image. 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: u1=0 on the left and the right edges of length h, with tractions T1=T2=0 on the top and the bottom surfaces of length l and T2=0 on the left and the right edges (figure 10a). (ii) Biaxial isometric tension: u1=0 on the left and the right edges of length h and u2=0 on the top and the bottom surfaces of length l. The tractions, T2=0 and T1=0, on the vertical and horizontal edges, respectively (figure 10b). Under both loading conditions, the fields (Embedded Image) are spatially uniform with no rotation of the material. All results are thus independent of the dimensions of the cell.

Figure 10

Schematic of the isometric (a) uniaxial tension and (b) biaxial tension boundary value problems under consideration. The coordinate system employed is also shown.

For uniaxial isometric loading, the time evolution of the normalized stress Σ11/σmax (figure 11a) and the corresponding evolution of the logarithmic transverse strain Embedded Image (figure 11b) 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 12a 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 12b). 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 12a at t/θ=0.5. However, fibres at ϕ≈90° are unconstrained and shorten at strain rate Embedded Image, 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 x2-direction. But since a small stress σ22 is generated in x2, 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 x2-direction. Increasing Embedded Image elevates the isomeric tension and decreases the transverse straining (figures 11 and 12). Moreover, with increasing Embedded Image 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 x2-direction. The effect of Embedded Image on the structural anisotropy is illustrated in figure 12, revealing that increasing Embedded Image causes fibres over a wider range of ϕ to be recruited into the tension generating apparatus, thereupon reducing the structural anisotropy.

Figure 11

(a) Time evolution of the normalized stress Embedded Image for uniaxial and biaxial isometric tension. (b) The corresponding time evolution of the logarithmic transverse strain, Embedded Image. The uniaxial isometric tension predictions are shown for three selected values of the normalized passive Young's modulus, Embedded Image. Unless otherwise specified, the material properties are taken as their reference values.

Figure 12

(a) Distribution of the stress fibre activation at selected times Embedded Image, as a function of the fibre angle ϕ, calculated for uniaxial isometric tension, with Embedded Image. (b) The steady-state distribution of the stress fibre activation as a function of the fibre angle, ϕ, for uniaxial isometric tension, with three selected values of Embedded Image.

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 11a) 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.

Figure 13

Schematic of the boundary value problem analysed to simulate the contraction of a square cell on an array of four posts. Note that the springs exert a force whose infinitesimal increments with respect to time are normal to the deformed perimeter of the cell. (a) Contraction of the cell constrained by the four supports and (b) contraction of the cell constrained by the four supports with simultaneous application of a localized tension on a patch on the right surface of the cell.

(a) Geometry and boundary conditions

A square cell with reference properties, side L=50 μm (thickness b=1 μm), is supported over a length Ls=5 μm at the four edges by identical elastic foundations of stiffness kE=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 Embedded Image on a cell surface with unit outward normal Embedded Image (in the current configuration) is given byEmbedded Image(5.1)where Embedded Image is the displacement rate along Embedded Image and Embedded Image is the magnitude of the traction vectorEmbedded Image. 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 Embedded Image and Embedded Image. In all the examples presented here, Embedded Image.

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

  1. The average stress fibre activation over all orientations,Embedded Image(5.3)

  2. The maximum principal value of the active stresses, σij, and the associated principal direction, measured as the orientation ϕp with respect to the x1-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 13a.

  1. A single activation signal with a decay constant of θ=1440 s applied to the cell at time Embedded Image (i.e. Embedded Image with all other cell properties fixed at their reference values).

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

  3. Four activation signals (each with a decay constant of θ=360 s) applied to the cell at times Embedded Image, 15, 30 and 45. Thus, Embedded Image 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, Embedded Image, and forces, Embedded Image (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 Embedded Image 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, Embedded Image. The long activation signal results in a higher concentration of stress fibres at Embedded Image, while at steady state (Embedded Image), 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.

Figure 14

Time evolution of (a) the average support displacement Embedded Image and (b) the support force Embedded Image for the cell sketched in figure 13a subjected to multiple activations and a single prolonged activation.

Figure 15

Distribution of the average stress fibre activation level Embedded Image for two of the cases considered in figure 14. (a) Two activations and (b) a single prolonged activation. The distributions of Embedded Image are shown at times Embedded Image. The distributions of the orientations, ϕp, of the maximum principal stress, Embedded Image, are also included as line segments (with the length scaled by the magnitude of the normalized stress Embedded Image). The solid circles show the original positions of the cell corners.

(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 13b 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 Embedded Image is applied in the x1-direction over a centrally situated patch, length Lf=5 μm, on the right surface of the cell (figure 13b) over duration, Embedded Image (where Embedded Image). Subsequently, the patch is fixed, with displacement rate Embedded Image prescribed in the x1-direction, over Lf.

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

The normalized force Embedded Image required to continue displacing the patch at a rate Embedded Image increases with time (figure 16). When the displacement rate is set to zero (at Embedded Image), the force begins to decrease before reaching steady state, Embedded Image, at time Embedded Image. The distributions of the average stress fibre activation levels at Embedded Image (figure 17a,b, respectively) with the orientations of Embedded Image 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 (Embedded Image), a high concentration of stress fibres continues to persist near the supports and at the site of the applied tension.

Figure 16

Time evolution of the normalized applied force, Embedded Image, for the problem sketched in figure 13b. The displacement rate Embedded Image for Embedded Image and the force is at its peak value at Embedded Image.

Figure 17

Distribution of the average stress fibre activation level, Embedded Image, for the problem sketched in figure 13b. The distributions are shown at times (a) Embedded Image and (b) at steady-state corresponding to Embedded Image. The distributions of the orientations, ϕp, of the maximum principal stress, Embedded Image, are also included as line segments (with the length scaled by the magnitude of the normalized stress Embedded Image). The solid circles show the original positions of the cell corners.

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.

  1. 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.

  2. 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.

References

View Abstract