Mechanosensitivity of the 2nd Kind: TGF-β Mechanism of Cell Sensing the Substrate Stiffness

Cells can sense forces applied to them, but also the stiffness of their environment. These are two different phenomena, and here we investigate the mechanosensitivity of the 2nd kind: how the cell can measure an elastic modulus at a single point of adhesion—and how the cell can receive and interpret the chemical signal released from the sensor. Our model uses the example of large latent complex of TGF-β as a sensor. Stochastic theory gives the rate of breaking of latent complex, which initiates the signaling feedback loop after the active TGF-β release and leads to a change of cell phenotype driven by the α-smooth muscle actin. We investigate the dynamic and steady-state behaviors of the model, comparing them with experiments. In particular, we analyse the timescale of approach to the steady state, the stability of the non-linear dynamical system, and how the steady-state concentrations of the key markers vary depending on the elasticity of the substrate. We discover a crossover region for values of substrate elasticity closely corresponding to that of the fibroblast to myofibroblast transition. We suggest that the cell could actively vary the parameters of its dynamic feedback loop to ‘choose’ the position of the transition region and the range of substrate elasticity that it can detect. In this way, the theory offers the unifying mechanism for a variety of phenomena, such as the myofibroblast conversion in fibrosis of wounds and lungs and smooth muscle cell dysfunction in cardiac disease.


Introduction
How cells sense is important to life: homeostasis necessitates sensors. The relationship between cell morphology and chemical properties of its environment have long been understood and documented. A plethora of sensory systems embedded in the cell membrane are able to detect specific signaling molecules, giving information to the cell to which it can react. Over the previous two decades there has been growing evidence to suggest that cells use these chemical pathways to sense the mechanical properties of their environment.
Mechanosensing is a process where at least one chemical reaction in the cell changes in response to a change in its mechanical environment. We specifically wish to distinguish two very different cases: when external forces are applied to the cell, and when the cell needs to measure the stiffness of its passive environment. We call the former case the 'mechanosensitivity of the 1 st kind'. There is a lot of research on this problem: the cell response to hydrostatic pressure [1,2], substrate topography [3] are just some examples; practically, there are implications for biotechnology [4] and cancer research [5,6]. Ultimately this kind of mechanosensitivity is well understood: there are many molecular processes that respond to an applied force and several different ones are employed by cells in different situations. It has been observed (for example) that mechanotransductive pathways can link compressive forces to translational and post-translational events by cytoskeletal interactions [7], that tyrosine kinases are part of the pathway of the cell response to shear stress [8], and that mechanical forces can alter ion channel permeability [9]. In each of these cases we find a remarkable protein design to transduce a mechanical stimulus into a chemical pathway that the cell can understand.
It is much more challenging to understand the mechanism of 'mechanosensitivity of the 2 nd kind', when the cell responds to substrate stiffness [10,11]. In order to detect a linear response function (elastic constant or viscosity), the cell has to apply a force of its own, and two separate measurements have to be carried out: of this force or stress as it is transmitted to the substrate; and of the displacement or strain in the deforming substrate. The "stiffness" is a ratio relating these two physical parameters. For this kind of measurement, the apparatus requires two points of application: to measure displacement a fixed reference point is always needed-to measure force one has to ensure a point of reaction as well. We need two fingers to squeeze a test object from two sides to determine how stiff it is-in the same way, all engineered devices ultimately have two points of action on the test sample. A mechanical sensor like the latent transforming growth factor-β complex we study in this paper (large latent complex, LLC) is effectively a single localised molecular protuberance outside the cell surface; it only has one point of application through which a pulling force is transmitted from the cell interior, and no a priori information about the reference point in the substrate [12]. An attempt to use a dynamical regime to overcome the need for a separate reference point in space also does not work: it is widely known that cells exert a constant force on their environment, not a dynamically varying force [13][14][15]. A recent paper [16] have formulated a new model, using the viscoelastic response of the substrate coupled with the stochastic process of rupturing the LLC by a constant pulling force generated by the cell. The result of this model is the predicted rate of free TGF-β release, expressed as a compact function of the pulling force and the substrate stiffness.
Cells in the body and in-vitro are normally supported biochemically and mechanically by the extracellular matrix (ECM), or its fragments (such as fibronectin) bound to artificial substrates. Experiments have observed that the stiffness of a substrate can stimulate differentiation of a fibroblast into a myofibroblast, or determine the future lineage of stem cells [17]: whether they will differentiate to bone, muscle, or nerve cells for example. How does the cell use proteins to sense stiffness and how does it convert a measurement it might be able to make into a chemical signal that the cell can understand and respond to? And what might a robust system which is designed to answer these two questions look like?
In this work we develop a coupled set of kinetic equations that describe a dynamical signalling feedback loop capable of interpreting the measurement of the substrate stiffness and convert it into the morphological changes the cell needs in reaction to such a substrate. We follow the earlier work by using the LLC family as a particular biological sensor, shown in [18] to be associated with the fibroblast to myofibroblast differentiation, and in [19,20] with the development of smooth muscle cells from embryonic stem cells, as well as their transformation following arterial injury or significant deposition of mechanically rigid plaque. We combine the result on the stiffness-dependent rate of latent complex rupture with the analysis of the coupled biological system to derive a set of equations to describe the dynamics of the protein pathway which forms the mechanosensing signalling loop. We compare the results and predictions to the stiffness-dependent behavior observed in cell-adhesion experiments, as well as in-vivo observations, and find a strong correlation of the model behavior and real systems.
In the discussion of the biological system and comparison to experimental results we will appeal most significantly to the fibroblast system for guidance. This is because a substantial volume of research has been carried out into the behavior of this system and because the model developed in [16], which we incorporate into our model, bases its arguments on the fibroblast system too. Besides, fibroblasts are the most prevalent cells in connective tissue and are vital in tissue repair within the body. When a break in the tissue occurs fibroblasts differentiate to myofibroblasts which act to pull the wound back together. The highly contractile myofibroblast form is produced by high expression of α-smooth muscle actin (α-SMA) fibres [21], the production of which has been shown to be induced by TGF-β [18,22].
The mechanosensing system is summarised diagrammatically in Fig 1. The latent complex of TGF-β is the stretch-dependent mechanosensor; it is representative of a larger family of mechanosensor complexes. There are three isoforms of the signalling molecule: TGF-β1, 2 and 3 with distinct effects in vivo, but their complexes are all produced in the same way: tightly but noncovalently bound as a dimer to their latency associated peptide (LAP-1,2,3) to form the 'small latent complex'. This forms a larger complex (which we are referring to as LLC) by disulphide links between certain domains in LAP and the latent TGF-β binding protein (LTBP), which in turn can bind to the ECM to anchor the complex on the substrate [23]. The large latent complex as a whole is secreted out of the cell, where it may connect LTBP to ECM components. Though no unambiguous binding partner for the complex has been found in the ECM [24], there is evidence in favour of strong covalent binding, likely to 8-Cys domains on LTBP [25] and it has been known to bind to certain ECM proteins such as fibrillin and fibronectin [26]. We shall assume that this bond is strong and not breaking during the processes we study.  [16]. The large latent complex of TGF-β binds to ECM proteins attached to a substrate, which may be deformable. The other end of LLC is attached to the cytoskeleton via an integrin complex that transmits the pulling force. If released, the free active TGF-β can bind to a set of receptors on the cell surface to initiate the signalling loop discussed in this paper.
On the other side, the LLC adheres to integrins on the cell surface. which are bound to the actin fibres making up the contractile cytoskeleton, allowing the pulling force generated by myosin motors to be transmitted to the complex, and through it-to the substrate. Upon mechanical rupturing of the complex, the active TGF-β is released. Once the LLC is spent, it can no longer function as the sensor and may or may not be replaced by a new latent complex.
Once TGF-β is released, it is free to diffuse. It could therefore bind to specific receptors on the cell surface [27]. Binding initiates the signalling loop in the cell. It may stimulate production of new LLC units to replace those that are spent, and also add more α-SMA on the inside of the cell [18]. α-SMA is produced to form more actin filaments and thus increase the pulling force F applied to the latent complex adhesion points. We shall assume a simple linear relationship between the concentration of α-SMA and the force F applied to LLC.

Rate of latent complex breaking
To describe the underlying physics of the system we start with the rupturing of LLC by an external pulling force F. On the far side this complex is attached to a viscoelastic substrate which we model by a simple Voigt model with an effective stiffness κ, that is directly related to the Young modulus E of the elastic medium via the classical continuum-mechanics relation: k % 4 3 pEx [28]. Here ξ is the characteristic length scale of short-distance cutoff in the Green's function of linear elasticity. We shall find that our results best match experimental observations when this length scale is taken as ξ % 100nm, which happens to be a characteristic mesh size of the extracellular matrix of crosslinked filaments. Then a rigid glass substrate, with a Young modulus E % 50 GPa, gives the maximum value of our κ % 21,000 N/m 2 . A typical muscle, organ or connective tissue (as well as synthetic gels) have E % 10 − 50 kPa, giving the range of κ % 0.004 − 0.02 N/m 2 . At the other end of rigidity scale, the brain tissue is very soft: E % 100 Pa, leading to the lower bound of κ % 4 Á 10 −5 N/m 2 .
Overdamped thermal fluctuations dominate the response of this series of mechanical elements. In equilibrium (i.e. with no force applied) the probability (rate) of spontaneous breaking of the latent complex is proportional to exp[−ΔG/k B T] with ΔG the bonding free energy of the large latent complex and k B T the thermal energy of the reservoir. When the pulling force is applied to the complex, it is also transmitted to the deformable substrate and the effective potential holding the latent complex together changes.
The rate of breaking of the latent complex attached to a viscoelastic substrate, under a constant pulling force, is given by the expression below (which corrects a typo in the original work [16]): where the shorthand non-dimensional notations used are: Here u 0 is the characteristic distance the latent complex needs to be stretched for rupturing. We could estimate the bond free energy ΔG holding the large latent complex together as for a few hydrogen bonds [29], ΔG % 10-16 kcal/mol, and the stretching distance before rupturing as the size of an aminoacid residue, u 0 % 0.3 nm. This is consistent with descriptions of force induced TGF-β activation by binding of α v β 6 to the latent complex [12]. In this case the range of the scaled (non-dimensional)k is between 4 Á 10 4 at the upper bound of rigid glass, 0.02 for a typical gel or muscle tissue, down to 8 Á 10 −5 at the lowest bound of E = 100 Pa. The effective diffusion constantD, which enters as a common factor in the expression for the rate of mechanical breaking k m , is the ratio k B T= ffiffiffiffiffiffiffiffi g 1 g 2 p reflecting the geometric mean of the damping coefficients (loss factors) of the substrate and the latent complex (see [16] for detail). We will find it necessary to estimate the value of this effective diffusion constant, which is determined by the effective damping coefficient of the viscoelastic substrate medium γ % κτ R (that is, the characteristic time of stress relaxation in the medium). There is a broad range of stress relaxation times found in different material systems, but for the sake of simplicity we shall take the characteristic Rouse time for a polymer network: τ R % 10 −5 s; whether it is a synthetic polymer gel or an organic tissue-this may serve as a first estimate (we must be aware that such a simple argument is likely to fail in the limit of very rigid glass and also in very soft gels/tissues with extra-long polymer strands). Then for a mid-range substrate stiffness of a typical gel or muscle tissue (E = 20 kPa), at room temperature we will haveD % 5 Á 10 À14 m 2 =s, a sensible value lower than the typical constant of thermal diffusion in water. We take the form ofD to be k B T/κτ R in the numerical analysis below.
In spite of its apparent complexity, the expression for the rate of mechanical breaking of the latent complex has only a few straightforward features: it describes the latent complex as having the physical bond energy ΔG in equilibrium, which is decreased when a pulling force is applied to it, making the basis for the Bell formula for the rate [30]: when the energy barrier is large, the applied force is small, and the substrate stiffness is low-Eq (1) does reduce to: However, we cannot use this Bell formula since we need to consider high pulling forces and moderate to low barriers that occur in practice. The correct expression Eq (1) has the appropriate behavior in the limits of large force, or small barrier, via the factor of inverse Δ F . On the other hand, Eq (1) also contains the effect of the deformed substrate, which is also subject to the same force F, transmitted through the latent complex while intact. This produces the 'enzyme effect' discussed in [16]: on increasing the pulling force and deforming the substrate into a new (deeper) elastic energy miniumum, the fluctuations of the point of adhesion become confined; this is expressed by the factor exp[−F 2 /2κk B T] in Eq (1). The other effect of the substrate viscoelastic response is to alter the non-exponential prefactor of the rate expression by adjusting the statistical weight of the intact latent complex in the full ensemble. The dependence of the breaking rate, k m , on the applied force for several values of substrate stiffness is illustrated in Fig 2. At sufficiently high substrate stiffness, expressed by the non-dimensional scaled parameter k given below the Eq (1), there is a maximum in the rate of breaking of the latent complex. So, initially, the positive feedback is expected-an increase in the pulling force increases the rate of latent complex breaking, leading to a further increase in force, etc. until the cell finds the equilibrium. However, at very low values ofk the maximum is at zero applied force. That is, on such soft substrates a pulling force actually reduces the rate of latent complex breaking causing the negative feedback. The outcome depends on the details of interactions between the key parameters in this feedback loop.

Protein interactions
Upon rupturing of one latent complex, an active TGF-β unit is released and becomes free to diffuse towards a specific receptor type on the cell membrane, and bind to it. The binding of TGF-β stimulates production of α-SMA [31] by a signalling pathway that we do not go into details of, except to know the ratio of how many α-SMA monomers are produced to each TGFβ bound. Given the linear relationship between α-SMA and force F applied to LLC, this stimulated production results in an increased force on the latent complex. The increase in applied force changes the rate of rupture of the complex, with the feedback nature depending on the regime of k m (F, κ).
We now seek to mathematically model the system of protein interactions described above. The three key proteins are α-SMA (labelled α), TGF-β (labelled β) and LLC (labelled λ), see To be specific with dimensionality, we consider the evolution of the number per cell of each of these proteins, α(t), β(t) and λ(t). The number per cell can be re-dimensionalised into a concentration by considering the typical mass of the protein in question and the volume of the cell-e.g. a typical fibroblast. The plot highlights the dramatic increase of the rate on stiffer substrates, and also the highly non-monotonic variation with force. The inset shows that on very soft substrates the rate of rupture is very low and in fact decreases further with the force. Note that the rate always diverges at f c ¼ 3=2, which is the maximum force that the latent complex can stand. Within our model there are certain assumptions about the state of proteins that we count in our variables. λ(t) represents only those LLC that are attached to both the cell and the substrate and have not ruptured. β(t) includes only TGF-β molecules that are released from their latent complex and are close to the cell surface, but not yet bound to a receptor. We can describe the evolution of these three variables in the following set of differential equations: Here the coefficients A λ and A α are, respectively, the number of fresh LLC complexes and α-SMA monomers produced for every instance of TGF-β binding to a receptor. Implicit within these two quantities are all the underlying protein pathways whose starting point is a TGF-β binding to its receptor, and whose end point is the production of a new latent complex and α-SMA subunit, respectively. k β is the rate of binding of TGF-β to the cell receptors. The key element of this model is the rate of the rupture of LLC, given in Eq (1) and discussed above.
We have to assume a linear relationship between α(t) and the pulling force F to be of the form: where z is the fraction of α-SMA monomers polymerised into actin filaments that are involved in the application of force on λ (this is expected to be a small number). f m is the force exerted by a myosin motor, which is a known parameter for each individual motor (* 2 pN), but may be a function of ATP availability in the vicinity of the adhesion complex [32]; n m is the number of myosin motors that act along each filament on average [33], a parameter that the cell can adjust within a wide range.

The dissipative sinks
We must also consider certain dissipative terms that affect all three of our dynamic variables, so that the system can find a homeostatic state; three dissipative terms are relevant, one for each of the three proteins. The λ sink is the saturation of integrin binding sites. We expect there to be only a certain number of sites over the surface of cell in contact with the substrate, at which the latent complex can attach to both the cytoskeleton and the substrate. This number may be large, but it is nevertheless limited, and this restricts the maximum value that λ can take. We consider a probability of a LLC attaching to the substrate and cytoskeleton which has the form: where N is the total number of available binding sites on the cell surface [34] and P 0 is the probability of a latent complex attaching to the substrate and cytoskeleton when N = 1 (i.e. there is no limitation on λ). This probability should be multiplied by the rate of production of LLC (A λ k β β) to give a total rate of production of λ, to replace the second term in Eq (3): Given the context of the problem, we expect the free TGF-β molecules to not only diffuse towards the cell receptors but also could drift away from the cell entirely. Therefore, a β sink due to diffusion away from the cell surface should be included. We consider this diffusion loss by treating TFG-β as a free Brownian particle diffusing away from a sphere (representing the cell surface). By considering the Gaussian probability for a molecule to be at a position r from the surface after time t, then setting r = 0 we find the concentration near the surface to decrease as P(r = 0) = ℓ(2πDt) −1/2 , where ℓ is some normalisation constant with units of length.
Our variable β(t) is then considered to be the probability of finding the TGF-β particle at r = 0. Initially we expect TGF-β to be present only at the cell surface, where it was released from the latent compolex. To normalise to unity at t = 0 we need to identify the radial lengthscale on which TGF-β exists. Because the probability tends to a delta function as t tends to 0, we assume that the radial lengthscale for normalisation is the diameter of TGF-β molecule: ℓ = 2R β .
Considering the probability as a volume fraction, we can write P(0, t = 0) = βV β /V tot where V β is the volume of one TGF-β subunit and V tot is the total volume occupied by all β at t = 0: a spherical shell at the radius of the cell (R cell ) of thickness 2R β . Therefore we want the contribution to the differential Eq (4) to yield: ffiffiffiffiffiffiffiffiffiffi ffi 2pDt p Þ, when no other terms are present. Differentiating gives a rate of diffusive loss: For our calculations it will be important to estimate the coefficient in this relation, which involves several known length scales; the values for D [35,36], R cell [37] and R β [38] are given in Table 1 below.
In the cell we always find that α-SMA decays over time, being used in various other processes. The rate of this spontaneous decay has been measured experimentally and also listed in the Table 1; we denote it as k α [31]: Combining all of the above together we arrive at the final set of three rate equations, which are now have balanced sources and sinks for the evolution of all three of our markers: Fixing parameter values The non-linear kinetic model of the signaling feedback loop has many parameters. Where available, values of experimentally derived parameters were always used in preference. The Table 1 gives the values of parameters that indeed have been determined experimentally. The value of k α was obtained from the line in Fig 5 of [31] representing decay of α-SMA without TGF-β in a cell anchored to a gel. By assuming exponential decay between day 1 and day 5, a value was approximated which fits this curve. The value of D is the diffusion constant for Green Fluorescent Protein (GFP), which is comparable in mass to TGF-β, and so it is expected that the values of both diffusion parameters are very close. No value of R β could not be found in the literature so it was approximated by assuming that GFP (a cylindrical protein) and TGF-β have similar volumes and that TGF-β is spherical. The diameter and length of GFP are 2.4nm and 4.2nm respectively. The values ofD, ΔG and u 0 for the large latent complex mechanosensor itself have been discussed earlier in the text.
Other parameters were not experimentally defined, so values that seemed physically reasonable were chosen and the dependence of system behavior on these values checked. These parameters are listed in Table 2.
The value for P 0 A λ was selected by assuming each TGF-β binding to a receptor stimulated production of approximately 1 or 2 new LLC. We expect all (or almost all) the LLC produced to be used in this mechanosensing system. We also assume that P 0 is close to 1. A value of 1.7 was selected to approximately reflect these assumptions. A α was selected to roughly reflect the high concentration of α-SMA in the cell. In order for the cell to react to a stimulus it will need to produce a significant amount of α-SMA monomers.
The value of z, the fraction of newly produced α-SMA that ends up polymerised into Factin and contribute to the pulling force on the adhesion complex, was chosen by considering the ratio of the volume of an α-SMA filament monomer (modelled as a sphere) to that of an α-SMA filament subunit (modelled as a cylinder). We define the length of a subunit as the step size of myosin which is 8 nm [39]. The diameter of actin filaments is about 8 nm [33]. The average diameter of an α-SMA monomer is about 4.8 nm [40]. This gives a ratio of 0.14. We then multiply this ratio by the fraction of α-SMA being used in this process to total α-SMA in the cell. This was estimated roughly to be 1/1500, but whatever this value-it is reassuring that this number is small since one does expect that only a small fraction of new α-SMA monomers would contribute to the increase of the pulling force on the point of adhesion.
The most difficult parameter for us to estimate is the rate of binding of free TGF-β to its receptors on the outer cell surface, k β . It is not experimentally determined and is difficult to define through arguments of plausibility. We fixed its value by requiring that the transition observed in steady-state values of λ, β and α (discussed in the following sections) corresponds to experimental results, so in a way k β is the only true fitting parameter of our model. In culture the fibroblast transition occurs at a Young modulus value E % 11 kPa [18]. Using k ¼ 4 3 pEx (as above, with ξ % 100nm) we find the corresponding dimensionless (scaled) elastic modulus: k transition % 5:9 Á 10 À3 . To reproduce this value in our results, the rate k β needs to be in the range shown in the Table 3.
Very recently, the binding rates of various TGF-β variants to cell receptors have been measured, returning rate constants in the range k on = 0.7-1.4 Á 10 6 M −1 s −1 [27]. The time constant of individual reaction that we use as k β is obtained by multiplying l on by the concentration of the TGF-β. There are several experimental measurements of this physiological concentration, giving the values in the range as 0.1-15 ng/ml [41][42][43]. For the sake of basic validation, let us take [c] = 1 ng/ml = 42 pM and obtain the product: k on [c] = 4.2 Á 10 −5 s −1 . Clealry the jigsaw pieces of this model fit together well.

Method of solving
All studies of the differential equations were done in Matlab using the 'ode23s' ordinary differential equations numerical solver package provided. Initial conditions were always set randomly, unless otherwise stated. For a given set of initial conditions the system always relaxed to Cell Mechanosensing the same steady state values. Fork ! 1 the solver would often take a very long time to produce a result, depending on initial conditions, because of the extreme stiffness of the system. Combinations of initial conditions andk that took a long time to solve were avoided for the majority of the calculation. However, the behavior was checked with one or two 'awkward' combinations to ensure that it agreed with other results and trends given in this paper. Note that to ensure λ reached steady state, a very long time span had to be used.

Evolution of variables
The variables α, β and λ evolve in time from various initial conditions to the final steady state. Evolution maps of the variables are shown in Fig 4 for lowk and a reasonably highk. The value for which we considerk to be 'high' or 'low' will be discussed later. It is clear that the dynamical system is robust and has a single stable attractor for a wide range of parameters (we have tested a large number of parameter values, not just the ones given in the Tables above).
The position of this fixed point attractor is dependent onk, that is, a single fixed point for each value ofk, and therefore provides a vehicle for mechanosensing of the 2 nd kind. Fig 5 shows the dependence of the time taken for β to reach steady state (relaxation time) oñ k. On the whole, on soft substrates the relaxation time varies from about 10 4.5 s to a little over 10 5 s (8 to 30 hours). Experiments have shown that cells react to TGF-β on time scales as long as 48 hours, the upper bound [41,44]-agreeing well with our results. It is interesting to observe how the time taken to reach the equilibrium steady state at highk increases substantially. At maximum possible values of substrate stiffness (k > 100) the average value is about 10 5.5 s (or 3-4 days). It is possible that this high values ofk are not practically achievable at all since rigid materials with a Young modulus E % 1-10 GPa have, in fact, a much smaller internal mesh size ξ. However, the overall trend of increasing the time to reach the homeostatic equilibrium with the substrate stiffness is a definite prediction of our theory.

Behavior in the steady state
The steady state dependence of each variable onk is of the highest importance in comparing our model to real biological systems and seeing how well it can describe reality. In a real cell, we would expect the concentration of TGF-β to determine the response of the cell. Different concentrations would result in different gene expression or different enzyme activation and overall a different behavior by the cell. We use the fibroblast system as the basis for discussion in this section. Fig 6 shows how the concentration of each variable behaves in the steady state. There are three distinct regions: the low and highk regions where all variables are unchanging, and a middlek region which is the transition region.
We first consider some peculiarities of the behavior of λ. In the highk region we see the value of λ tends to very close to zero, suggesting that as soon as LLC binds to the integrin and substrate-it ruptures, so that on average there is none (or very little) present. This is a strange result and is certainly not intuitive. The reason for this is that in the highk region focal adhesions form (illustrated in Fig 7) where many LLC cooperate together in a combined structure. This is not accounted for in our model because we assumed that λ represented only LLC attached individually to the substrate and integrin. Therefore we can deduce that only a small fraction of total LLC exist as individual elements in the highk limit-a deduction supported by experiment [45] showing that focal adhesion size increases on stiffer substrates. The true number of LLC may be much higher than that implied by the value of λ.
Examining the behavior of α and β in steady state, we see that the transition occurs in the region 10 À4 <k < 10 À1:5 . The corresponding Young's modulus where the fibroblast transition occurs lies in between these two extremes, as required by our fitting parameter.
Over this transition β increases from 6cell −1 to a final value of 108cell −1 . These values correspond to real concentrations of about 0.48ng/ml and 8.6ng/ml, respectively-taking TGF-β molecular mass as 25kDa [44] and assuming the fibroblast is a spherical cell with a radius of 5μm (a big approximation given the myofibroblast is strongly contracted). How does this compare to real systems? TGF-β concentrations as low as 0.1ng/ml were found to be sufficient to  induce a transition in endothelial cells to myofibroblasts with a maximal response at 1ng/ml [43]. Other reports in the literature [41,42] found stimulation of the fibroblast to myofibroblast transition or in dendritic cells for concentrations of TGF-β between 0.1ng/ml and 15ng/ ml.
The breaking rate k m ðf ;kÞ and applied force in equilibrium Finally, we examine the relation of force and rate of rupture of LLC. Qualitatively it would seem reasonable to expect the steady state to correspond to the maximum in k m (Fig 2) where the positive feedback changes to negative. Surprisingly we do not find this intuitive relation between system variables and k m . Fig 8 shows the relation between the steady state value of α (directly proportional to force) and the rate of rupture of LLC. As can be seen, system variables do not relax on the maximum in k m . Instead there are three distinct types of behavior. At very lowk we see the steady state force is beyond the maximum in k m which is atf % 0 (so the cell is not predicted to reduce its pulling to zero, as might have been assumed). At highk on rigid substrates, the steady state force occurs before the maximum, so in spite of its restructuring and formation of more stress fibers-the coupled signalling loop stops the cell short of reaching the maximum rate of latent complex breaking. At intermediate values, in the 'transition region' ofk, the steady state force rests somewhere after the peak of k m .
We find that the transition in λ from a high steady state concentration at lowk to a low concentration at highk occurs between the low and intermediate cases discussed in this section. No particular effect in λ is observed when the steady state value of α moves past the maximum in k m . As far as the system has been tested, it seems that the maximum in k m has very little part to play in the dynamics of our coupled system variables.
Using n m to adjust the transition region The system proposed is clearly not sensitive to the substrate stiffness in the low or highk regions. However, it has been observed that cells respond to elasticities out of the range our model would suggest is possible [17]. We found that changing n m (the 'number' of active myosin motors on the actin fiber) does not affect the steady state values of any variables away from the transition region. But it does affect the position of the transition region versusk-as seen in Fig 9. If we allow the cell the flexibility to vary n m as well as the applied force, then n m could be adjusted so that the transition region scans over many orders of magnitude ofk. Myosin motor dynamics can be regulated by de/phosphorylation in-vivo [46]. Therefore, as long as the cell has a method to record n m and one of the variables, then it could dynamically vary n m to find the elasticity of its substrate.

Conclusion
This work was a natural continuation of the earlier study [16] of the rate of rupture of LLC. Incorporating their model into a full dynamically coupled system allowed us to create a set of differential equations with the mechanical sensing system encoded in this rate. The behavior of the equation is entirely dependent on the applied force and the elasticity of the substrate,k, and we have checked by starting from a large variety of random initial conditions that the equilibrium (steady state, homeostasis) of the cell is unique and specific on each substrate. We found that the system always returns to the steady state independently of initial conditions (a single stable attractor) and there appears to be a certain characteristic time scale of relaxation towards this equilibrium, regardless what the cell parameters have started with, of the order of 10 4.5 s (i.e. 8-9 hours), which is a broad agreement with in-vitro observations of cell deposited on various substrates. The choice of model parameters may appear the most controversial aspect of our work (since there are so many of them and some are little-known), yet we have taken great care to fix most of these parameters at specific values found in experiment. It is clear that for any particular cell system (osteoblast, fibroblast, smooth muscle, glial-or indeed stem or various progenitors) the values of these parameters would be subtly different, and the cell response range vary accordingly. But the range of this variation is not great, and we are assured that the qualitative nature of the predicted response will remain the same.
To summarise the predictions of our theory: 1. We study a response of a 'model cell' when it is placed on a specific substrate, with given elastic modulus and loss factor, or in other words-stress relaxation time. The principal sensor is the large latent complex of TGF-β, which we use as a representative model. Following [16], we found the compact expression for the rate of latent complex rupture, k m (F, κ) in Eq (1), which is determined by the parameters of the substrate and the steady pulling force from the cell. This breaking rate has a complex non-monotonic behavior with changing parameters, and it is central to the whole process of mechanosensitivity of the 2 nd kind.
2. On release of the signal molecule (active TGF-β) after the irreversible breaking of the latent complex, it could bind to a specific receptor on the cell to initiate the feedback loop-or it could drift away from the cell. We have not yet followed an obvious conjecture, that these 'lost' TGF-β signal molecules may provide an element of quorum sensing in a confluent tissue or cell colony. The concentration of free TGF-β released from the broken latent complexes (our variable β) has specific values in homeostatic equilibrium.
3. On placing the cell from any previous environment onto the given substrate (which we modelled by scanning a broad range of initial conditions), the cell remodels itself until it achieves the 'equilibrium' morphology specific for this substrate. The most obvious and experimentally tested signature of this is the concentration β in the vicinity of the cell, which is low/vanishing on very soft substrates, and high constant on mechanically rigid substrates (which is what has been registered in many experiments on very different cells and environments). The transition region in between these two limits is where the cell sensitively responds to small changes in stiffness by altering the equilibrium values of its parameters.
4. The adjustment of cell morphology to its substrate stiffness happens via the expression of additional α-SMA (stimulated by the TGF-β binding), which is then incorporated into additional F-actin and stress fibers-and provides more strong pulling force F(α) on the latent complexes in points of adhesion. Again, in close agreement with observations, we find that actin production is low on soft substrates and high on rigid substrates, with a transition region in between. The fibrosis in tissues is therefore linked to the stiffness of the environment the fibroblasts or smooth muscle cells are in contact with.
5. Our final kinetic variable, the concentration of intact large TGF-β latent complexes (λ) shows an interesting an unexpected response to stiffness change. On soft substrates, with low α (so the cell does not pull too hard on its adhesion points) and low β (reflecting a low probability of LLC breaking)-there are a lot of adhesion points around the cell: the equilibrium level of λ is high. However, on rigid substrates when the cell pulls with a large force, we find λ dropping to a very low value. Note that since the level of β is high in this region, the breaking of latent complexes occurs very frequently-almost immediately after they have been replenished by the cell, so there are almost no active adhesion points in operation around the cell surface. This may sound contradictory, but we interpret it as the commonly observed state when just a few large focal adhesion complexes hold the cell on a rigid substrate-and we have now discovered that the cell surface outside of these few focal adhesions is free from the 'singular' adhesion points. This is in contrast with the same cell on a soft substrate when it has a homogeneous concentration of singular adhesion points (high-λ state). Our theory does not describe the collective effects that take place in large focal adhesions: certainly both the pulling force and the latent complex rupturing are very different there.
We have tried to make sure the stiffness transition region where steady state concentrations of α and β rise by many factors, and λ falls significantly, corresponds to observations for fibroblasts. Since most of the parameters we use are independently fixed, the critical role is played by the binding rate of TGF-β to its receptor k β , which we assume is not known. The value of k β % 5 Á 10 −5 s −1 was required to bring this transition region to the Young modulus E * 10-15 kPa required to observe fibroblast to myofibroblast transition. However, we have shown that the cells have several mechanisms to adjust this transition region even for exactly the same LLC adhesion parameters: for instance by varying the number and the efficiency of myosin motors (the parameter n m , Fig 9), or by channeling a different amount of newly produced α-SMA to the adhesion pulling filaments (the parameters z and k α ). Therefore we would expect that glial or neuron cells would 'move' their stiffness transition to a much lower Young modulus range, while osteoblasts may benefit from their region of primary sensitivity to lie at a higher Young modulus as appropriate to their target environment. Perhaps this choice is made at the stem level and is at the core of subsequent differentiation: by dynamic variation of n m by de/phosphorylation the cell can 'scan' the transition region over many orders of magnitude of elasticity.
We found that the system had no particular dependence on the prominent and characteristic peak of the LLC rupture rate k m (F, κ), defying our intuitive expectation for how the system might respond. It seems that the general pattern of behavior in homeostatic equilibrium is independent of many parameters-at least to the extent that it was tested. Further work may seek to answer some of many questions posed by this paper. Does it matter that the peak in k m seems to not have the importance we expected? How could we include the effect of focal adhesions on LLC and would it change the behavior of the model dramatically?