Coupling biochemistry and mechanics in cell adhesion: a model for inhomogeneous stress fiber contraction

Biochemistry and mechanics are closely coupled in cell adhesion. At sites of cell-matrix adhesion, mechanical force triggers signaling through the Rho-pathway, which leads to structural reinforcement and increased contractility in the actin cytoskeleton. The resulting force acts back to the sites of adhesion, resulting in a positive feedback loop for mature adhesion. Here we model this biochemical-mechanical feedback loop for the special case when the actin cytoskeleton is organized in stress fibers, which are contractile bundles of actin filaments. Activation of myosin II molecular motors through the Rho-pathway is described by a system of reaction-diffusion equations, which are coupled into a viscoelastic model for a contractile actin bundle. We find strong spatial gradients in the activation of contractility and in the corresponding deformation pattern of the stress fiber, in good agreement with experimental findings.

adhesion-related processes, cellular behaviour is not only controlled by biochemical cues, but also involves many physical determinants like the structural organization of the extracellular matrix and the cytoskeleton or force generation through molecular motors [5,6,7,8].
For example, it has been shown that the stiffness of the extracellular environment determines migration of tissue cells [9] and differentation of stem cells [10]. In particular, these celluar responses have been found to depend on the ability of the cells to contract their environment with actomyosin contractility and to convert this mechanical process into a biochemical signal. Although these processes are essential for such important situations like tissue functioning or cancer cell migration, theoretical models describing the coupling between biochemistry and mechanics in cell adhesion are still rare, albeit essential for a future systematic understanding of how multicellular organisms organize themselves.
Cell adhesion is closely related to the actin cytoskeleton, whose organization is central in determining the structural properties of cells. In cell culture with stiff substrates, the actin cytoskeleton tends to organize in stress fibers, which are bundles of actin filaments tensed by myosin II molecular motors [11]. Stress fibers usually end in focal adhesions, which are integrin-based adhesion contacts which can grow to a lateral size of several microns [12,13]. On their cytoplasmic side, focal adhesions recruit more than 90 components (mostly proteins) which physically reside in the adhesion structure [14]. In 1992, Ridley and Hall published a landmark paper demonstrating that the assembly of stress fibers and focal adhesions is regulated by a small GTPase called Rho [15]. Rho has many isoforms, but the one mainly associated with focal adhesions is RhoA, which for simplicity in the following we refered to as Rho. In a companion paper of the same year, Ridley and Hall showed together with coworkers that another small GTPase called Rac stimulates the formation of lamellipodia as they appear in cell migration [16]. The main isoform associated with focal adhesions is Rac1 which for simplicity in the following we will refer to as Rac. While Rho mainly acts through activation of actomyosin contractility, the main effect of Rac is activation of actin polymerization, in particular activation of the actin nucleation factor Arp2/3. It has been reported later that activation of Rac downregulates Rho, leading to disassembly of stress fibers and focal adhesions [17]. In many situations, Rho and Rac can be regarded as antagonists, switching the cytoskeleton between different structural states [18]. They are part of a larger family of small GTPases, called the Rho-family, which for example also includes Cdc42, which stimulates the formation of filopodia and maintains Here we focus on the role of Rho as stabilizing factor for mature adhesion. During recent years, it has been shown that Rho is the central component of a biochemical-mechanical feedback loop which regulates mature adhesion. In detail, it has been shown that application of force on focal adhesions triggers their growth in a Rho-dependent manner [21] (reviewed in [13,18]). Two main downstream targets of Rho leading to stress fiber formation have been identified. The formin mDia leads to actin polymerization, while the Rho-associated kinase ROCK leads to phosphorylation of myosin light chain and thus to increased motor activity. Together these effects lead to formation of and contractility in stress fibers and therefore to increased force levels at focal adhesions. In this way, a positive feedback loop is closed for upregulation of mature adhesion characterized by focal adhesions and stress fibers.
This biochemical-mechanical feedback loop is schematically depicted in Fig. 1. An essential part of this feedback loop is the growth of focal adhesion under force, which recently has been the subject of different modelling approaches [22,23,24,25,26,27] (reviewed in [28]).
However, these models have focused mainly on the mechanical and thermodynamic aspects of the growth process, neglecting the interaction of mechanics and biochemical signaling. The positive feedback loop between contractility and growth of adhesions has been modelled before in the framework of kinetic equations, but without addressing the details of force generation and its regulation by signaling pathways [29]. Similar kinetic equations have been used to model the antagonistic roles of Rho and Rac in cell adhesion, but again without addressing the details of force generation and regulation [30]. Recently force generation has been addressed in more detail in a model for whole cell contractility and stress fiber formation [31,32]. However, no details of the signalling pathway have been modelled except for an unspecified activation signal.
Because the actin cytoskeleton is very dynamic and interacts with many different molecular factors, including actin-associated proteins and molecular motors, it is very difficult to model its mechanical properties in a general way. However, modelling becomes feasible if one focuses on one of the well-characterized states of the actin cytoskeleton, for example the lamellipodium or stress fibers. Because here we are mostly interested in mature cell adhesion in culture, we will focus on the latter case. Modelling stress fibers can be approached from different perspectives. An obvious starting point are their common characteristics with muscle fibers, which is a linear sequence of sarcomeres, each containing around 300 myosin II molecular motors working collectively together as they slide the actin filaments relatively to each other. This field has been pioneered by the Huxley-model [33], which later has been modified in many regards, e.g. in regard to filament extensibility [34] or by a detailed modelling of the myosin II hydrolysis cycle [35]. In contrast to muscle fibers, stress fibers are more disordered and a complete description therefore requires a model for their assembly process from polar filaments interacting through molecular motors. Such a description has been achieved in the framework of a phenomenological theory which however does not model the details of the underlying motor activity [36,37]. This theory does predict different dynamical states of the system, including a stationary state of isometric contraction as observed in stress fibers.
Although being less ordered than muscle on the level of electron microscopy, stress fibers do exhibit a periodic organization. Fig. 2 shows experimental data for fibroblast adhesion to a stiff substrate [38]. The image of a whole cell shown in Fig. 2A reveals a banding pattern in the stress fibers. The green regions correspond to the actin crosslinker α-actinin while the red regions correspond to myosin II molecular motors. Non-muscle myosin II is known to assemble into bipolar minifilaments consisting of 10-30 myosins, as depicted in the cartoon of a stress fiber in Fig. 1. In contrast to muscle, the striation pattern of stress fibers shows considerable variability along the length of a stress fiber. In particular, it has been observed that upon stimulation of contraction with the drug calyculin A, only the sarcomeres in the periphery (close to the focal adhesions, Fig. 2B) shorten, while those in the center (close to the cell body, Fig. 2C) elongate. The exact time course of the striation width is shown in Fig. 2D and E for control and stimulation experiments, respectively. This demonstrates that activation of contractility leads to strong spatial gradients in the striation pattern. Thus it is not sufficient to model only one sarcomeric unit as it is typically done for muscle [33,34,35].
Rather at least a one-dimensional chain of such sarcomeres has to be considered.
Here we introduce a new one-dimensional model for stress fibers which essentially models them as linear sequence of viscoelastic Kelvin-Voigt bodies, whose stationary state is determined by the elastic part. The action of the molecular motors inside each sarcomeric unit are included on the level of a linearized force-velocity relationship. The signaling pathway is modeled as a system of reaction-diffusion equations. For this study, we have conducted an extensive survey of the relevant literature and have collected the measured rate and diffusion constants in such a way that they now can be used for mathematical modelling.
The coupling between the biochemical signaling pathway and the mechanical stress fiber model proceeds by introducing a spatially varying fraction of active molecular motors. The local activation level is thereby determined by the outcome of the signaling pathway. A continuum limit of the mechanical model for many sarcomeric units in series results in a partial differential equation with mixed derivatives. The whole system of reaction-diffusion equations for the signal transduction and the partial differential equation for the mechanical part can be solved simultaneously. By feeding the force resulting from the mechanical model back into the activation of the signaling pathway, we obtain for the first time a model for the closed biochemical-mechanical feedback cycle described above. As a first application of our model, here we show that it predicts heterogenous contraction of stress fibers, in good agreement with experiments. In general, our work shows how models for the coupling of biochemistry and mechanics can be deviced in a meaningful way.
The paper is organized as follows. In Sec. II we start with our model for the part of the Rho-pathway which is relevant for our purposes. We then introduce our mechanical model for one sarcomeric unit of the stress fibers in Sec. III and its continuum limit for a whole stress fiber in Sec. IV. In Sec. V our model predictions for inhomogeneous contraction upon activation of contractility by calyculin are compared to the experimental results. Finally we conclude in Sec. VI.
In this study we concentrate on stress fibers and their regulation by the Rho-pathway.
In Fig. 1 we introduce a coordinate system for our one-dimensional model: the stress fiber extends along the positive x-direction and the endpoints at x = 0 and x = L correspond to two focal adhesions. Because the two focal adhesions are treated as equivalent, our model has inflection symmetry around x = L/2. In Fig. 3, we schematically depict the biochemical part of our model in the spatial context of the focal adhesion at x = 0 (by symmetry, the same description applies to the one at x = L). Three compartments have to be considered: the focal adhesion, the cytoplasm and the stress fiber. In our model, each of these compartments corresponds to one or two important biochemical components. The reaction pathway is a linear sequence of activating or inihibitory enzyme reactions initiated at focal adhesion, transmitted through the cytoplasm by diffusion and resulting in spatially dependent myosin activation in the stress fiber. In the following we discuss each reaction step in detail and show how these processes are translated into reaction-diffusion equations.
The abbrevations used for the biochemical components are compiled in Tab. I, together with a short description of their functions. The model equations are summarized in Tab. II.
Our modelling of the biochemical signaling pathway starts with the activation of Rho at focal adhesions. The Rho-protein has a lipophilic end that serves to anchor it to lipid membranes [39]. Complexation with guanine nucleotide dissociation inhibitors (GDIs) shields the hydrophobic parts of the Rho-protein and make it inactive as well as soluble in the cytoplasm [40]. It is expected that Rho is released from these complexes at focal adhesions.
More importantly, focal adhesions are known to recruit different Rho-GEFs, thus activating Rho at the focal adhesions. The active Rho is then able to bind Rho-associated kinase (ROCK) and thereby activate its kinase activity [41]. Since the active ROCK is bound to Rho-GTP, we assume in our model that these components are not diffusible but are localized to the focal adhesions. For the same reason, we neglect the direct interaction of ROCK and myosin which has been reported to occur in vitro. Active ROCK phosphorylates the diffusible myosin light chain phosphatase (MLCP) at its myosin-binding subunit (MBS).
MLCP and its antagonistic partner myosin light chain kinase (MLCK) are the main regulators of myosin contractility in the context we are interested in. Both enzymes interact with the regulatory myosin light chain subunit (MLC) of the myosins. Depending on the respective activities, MLC gets either phosphorylated or dephosphorylated. MLC, in turn, controls the myosin binding to actin filaments. Only if MLC is phosphorylated myosin is able to bind actin filaments and perform its ATPase cycle that converts chemical energy into mechanical work, e.g. contraction of stress fibers [42]. By phosphorylating MLCP, ROCK effectively enhances the phosphorylation level of MLC. In this way, Rho-activation can lead to an increase in myosin contractility in the stress fiber.
Our model for the biochemical reaction-diffusion system assumes that each enzyme stimulation follows Michaelis-Menten kinetics [43]. In Michaelis-Menten kinetics, production first increases linearly with educt concentration and then saturates at a maximal production velocity if educt concentration exceeds the value set by the Michaelis-Menten constant. The precise molecular process corresponding to the conversion of force into a biochemical signal at focal adhesions has not been identified yet. However, it is expected that mechanical forces exerted onto focal adhesions eventually initiate the loading of Rho-GTP leading to ROCK activation. For lack of information we therefore lump the focal adhesion associated processes into one equation that effectively describes the conversion of ROCK into its activated form (presumably complexed with Rho-GTP). The mechanical force F b that stimulates the activation is treated as enzyme in the framework of Michaelis-Menten kinetics: The variable ROCK denotes the activated form of ROCK and we assume that the overall concentration of ROCK is constant at ROCK tot . The force exerted by the stress fiber onto the focal adhesion, F b (t), stimulates the conversion of ROCK into its activated form with maximum velocity r 1 F b (t) and Michaelis-Menten constant K 1 . The parameter r 1 is equivalent to a rate constant but relates mechanical force to a chemical reaction. For this reason the units of r 1 are given as [nM/s nN ]. The force F b (t) will depend on the stress fiber deformation. The second term accounts for the degradation of activated ROCK to its inactive form, with maximum velocity V −1 and Michaelis-Menten constant K −1 . Since we expect ROCK in its active form to be associated with focal adhesions, we omit diffusive contributions to this equation.
One main effector of ROCK is MLCP which we regard as a diffusible compound leading to a reaction-diffusion equation: Here, the variables MLCP-P and MLCP denote the phosphorylated and unphosphorylated form of myosin light chain phosphatase, respectively. The first term accounts for the dephosphorylation of MLCP-P with maximum velocity V −2 and Michaelis-Menten constant The second term allows for the diffusion of the phosphatase with diffusion constant D. The phosphorylation level of MLCP is also regulated by the active form of ROCK which catalyzes the reverse reaction, that is the conversion of the phosphatase into its phosphorylated form. However, ROCK is only active in the vicinity of focal adhesions located at each end of the stress fiber. Therefore this source term can be incorporated into the boundary conditions for Eq. (2), in the sense that the diffusive flux into the boundary has to balance the conversion into its inactive form: The same relation, but with inverted sign is valid at the other end at x = L, compare the overview in Tab. II. The reaction is again modelled with Michaelis-Menten kinetics, where By allowing only the ratio of the phosphorylated fraction to vary, we assume that the overall amount of myosin in the stress fibers is fixed. MLC is phosphorylated by MLCK with a maximum velocity V 3 = r 3 MLCK and respective Michaelis-Menten constant K 3 . Here we assume that the concentration of MLCK is constant within the cell. The kinase is antagonized by MLCP that dephosphorylates MLC with a rate constant r −3 and Michaelis-Menten The factor I is an inhibition parameter defined below. Since MLCP has spatial dependent source terms and is diffusible, the inhibition of MLC by the phosphatase will vary in space.
To complete the biochemical modelling we have to specify how the induction of calyculin is treated in our model. Calyculin is an inhibitor of MLCP and thereby enhances the phosphorylation level of MLC. We model the interaction of calyculin with its target MLCP as a competitive inhibition leading to the additional factor I in the last term of Eq. (4) [44].

III. SARCOMERIC UNIT OF STRESS FIBER MODEL
A minimal model for stress fibers has to take into account not only the viscoelastic but also the contractile properties of the fiber due to myosin motor activity. For the mechanical response of a sarcomeric unit we take the usual Kelvin-Voigt model for viscoelastic material [46]. It consists of a dashpot with viscosity η and a spring of stiffness k connected in parallel.
These two modules represent viscous and elastic properties of the material, respectively.
The Kelvin-Voigt model is the simplest viscoelastic model which in the stationary state is determined by elasticity, in contrast to the Maxwell model, which flows in the stationary limit. Thus the Kelvin-Voigt model is the appropriate choice for stress fibers, which can carry load at constant deformation over a long time. In order to cope with the contractile behavior of stress fibers, we introduce an additional contractile element that represents the activity of motor proteins. For illustrative reasons we first derive the governing differential equation for such a viscoelastic and contractile element (vec-element) before we proceed to the stress fiber model. This ansatz is similar to the two-spring model which we have introduced before to explain the physical aspects of rigidity sensing on soft elastic substrates [47].
The considered vec-element is depicted in Fig. 5. The properties of the additional contractile element is given by the specific force-velocity relation of the molecular motor. For simplicity, we use the linearized relationship: F m is the actual force exerted by a motor moving with velocity v. v 0 is the zero-load velocity and F stall is the stall force of the motor, that is the maximal force allowing motor movement. In the following description for the stress fiber, the stall force will depend on the active fraction of the myosin, compare Eq. (13), such that myosin contractility may vary spatially. The sum of all internal forces, F, exerted by the vec-element reads The crucial point is that the motor force F m (v) is related to the contraction velocity v which The minus sign comes from the fact that the directed motor movement causes a relative sliding of the anti-parallel orientated actin filaments leading to contraction of the element with velocity v. Thus the displacement u(t) becomes negative upon motor activity. If there are no external forces, all internal forces have to balance, F = 0, and we arrive at the following differential equation for the displacement u(t): The term originating form the motor activity decomposes into two contributions.
Here τ = η e /k corresponds to the relaxation time of the Kelvin-Voigt model with an effective viscosity η e = η + F stall /v 0 . Thus the system appears to be more viscous for stronger and slower motors and subsequently the motor activity slows down any relaxation due to perturbations of the system. The steady state deformation u(t → ∞) = −F stall /k is determined by the ratio of stall force and spring stiffness.

IV. CONTINUUM VERSION OF STRESS FIBER MODEL
We model a stress fiber as a string of vec-elements whose viscoelastic and contractile properties may vary spatially. For example the contractile strength of the vec-elements will vary spatially because the biochemical signal will cause different myosin activation levels along the fiber. Starting from a discrete description depicted in Fig. 6 we derive the governing continuum equation. Similar to the preceding discussion, the force F n on the site n is the sum of spring forces, viscous drag and the forces built up by the motor proteins: Here, we allow that the stiffness of the spring, the viscosity as well as the motor force vary spatially. In order to deduce a continuum description we expand the functions η, k, u and F m at x = na assuming that these functions are smooth within the small distance a. Keeping only leading order terms yields: Note that the leading differential operator ∂ x acts on η and u. The same holds for the second term. If k does not vary spatially, it simplifies to k∂ 2 x . Like for a single vec-element, we argue that the motor force F m depends on the displacement u(x, t). The contraction ∆ n within the n-th element generated by the respective motor is given by ∆ n = −(u n − u n−1 ). The contraction velocity is therefore v(x, t) =∆(x, t) ≈ −a∂ t ∂ x u(x, t). The found expression for the velocity is inserted into the force velocity relation Eq. (5) leading to: In contrast to Eq. (5), here the stall force is not constant but depends on the phosphorylated fraction n(x, t) of MLC along the stress fiber. This comes from the fact that along a myosin minifilament and depending on MLC phosphorylation, a larger or smaller fraction of myosin heads is able to bind to actin and perform ATP-cycles. The more myosin heads are active the larger the maximum force that the bundle can exert to the actin filaments. In our model we regard the ensemble of myosins within a cross-section of a stress fiber as one large contractile unit with an effective stall force that depends linearly on the active fraction n of myosin heads: The effective stall force, F stall (x, t), would reach the maximum force F max if all myosins within this cross section would be working (n = 1). In the following we set F max = 50 nN which will result in boundary forces exerted by the fiber of about 5 nN which corresponds to typical values observed in experiments [45]. Eq. (11) together with Eq. (12) and Eq. (13) lead to the final model equation for the stress fiber: where the effective viscosity η e (x, t) = η(x) + F stall (x, t)/v 0 , similar to the findings in Eq. (8), although here the viscosity varies spatially. Interestingly, only the variation of the motor force appears on the right hand side of the equation. As a consequence, a homogeneous motor activity will not contribute to the displacements within the string.
To obtain some intuition for this equation, assume that the spatially varying stall force is given e.g. by the steady state solution n ss of the reaction diffusion system, depicted in Fig. 4, such that F stall (x) = F max n ss (x). For this simplifying case where the stall force does not vary in time, Eq. (14) can be integrated and the time dependent solution for u(x, t) is given by: Here, we have set τ (x) = η e (x)/k(x), the typical relaxation time with which perturbations decay at a certain position. E.g. the initial conditions ∂ x u(x, t 0 ) can be regarded as perturbations to the steady state and they decay with exp(−t/τ (x)). The three integration constants can be identified as the displacement at the left boundary, u(x 0 , t), the force exerted to the left boundary, F b (t), and the initial strain along the fiber, ∂ x u(x, t 0 ). They are determined by the boundary and initial conditions. Experiments by Peterson et al. [38] are arranged with cells on stiff substrates to which the ends of the stress fiber are connected by focal adhesions. Therefore, the appropriate boundary conditions are fixed ends for the fiber, namely u(x 0 , t) ≡ 0 and u(x e , t) ≡ 0. From the second condition one is able to calculate the missing integration constant F b (t) for any initial condition ∂ x u(x, t 0 ). The force on the left boundary F b (t) is given as solution to an inhomogeneous Volterra equation of the first kind: with the kernel and inhomogeneous part g(t) dependent on the initial condition ∂ x u(x, t 0 ): In order to solve the integral equation Eq. (16) for F b (t), we calculate its time derivative, leading to: This equation yields an explicit expression for the initial force F b at t=0: By inspection of the kernel Eq. (17) and the inhomogeneous part Eq. (18) one finds that the initial force onto the boundary has a finite value, even for the initial condition ∂ x u(x, t 0 ) = 0.
Eq. (19) also yields an iteration rule for the time course of F b (t), by applying a quadrature, where F b (t 0 ) from Eq. (20) is used as a starting value. The solution for F b (t) is shown in Fig. 7. The boundary force is rising from its initial value and then quickly saturates at about 4.8 nN. The result for F b (t) can then be set into the general solution (Eq. (15)) for the displacement u(x, t) along the fiber. Fig. 8 shows this solution, by using the steady state activation level for the myosins (F stall (x) = F max n ss (x)) shown in Fig. 4 and assuming the initial condition ∂ This relation now connects the biochemical signaling to the mechanical deformation of the stress fiber. Thus, the coupled system of reaction equations, Eq. (1) to Eq. (4), and the mechanical equation Eq. (14) have to be solved simultaneously. This can be done numerically by using the MATLAB algorithm "pdepe". The whole system of equations and the used parameter values are summarized in Tab. II and Tab. III, respectively.
By doing a steady state analysis we find, that this system of equations exhibit two stable steady states for the used parameter values: the first state is characterized by a generally low activation level n ss (x) of the myosin motors resulting in marginal boundary forces whereas in the second state myosin motors are non-uniformly activated and the exerted forces reach a few nN. This bistability is characteristic for a positive feedback system [3]. The first "nonactive" state would correspond to cells that failed to establish mechanical stress whereas the second "active" state correspond to cells that are well adhered to the substrate. In order to simulate the drug experiments by Peterson et al. [38], we start with the system residing in this "active" state and then, at t=0, we perturb the system by turning on the stimulation with calyculin. This is modelled by switching instantaneously the inhibition parameter from We have to stress that the presented stress fiber model is continuous, thus the model cannot distinguish between α-actinin bands or MLC bands. The color code in Fig. 10 and Fig. 13 is therefore arbitrary. We also derive the local relative change of density within the fiber which is given in general as the negative trace of the strain tensor δρ rel = (ρ − ρ 0 )/ρ 0 = −tr(u i,j ).
Since the model is one dimensional this simplifies to δρ rel = −∂ x u(x, t), plotted in Fig. 12.
The local relative change of band width at a certain position within the fiber, is then simply given by: t). The figure shows that the inhomogeneous motor activity causes a contraction of the bands up to about 55% close to the fiber ends (the relative change in density is positive), whereas the pattern expand up to 15% at the middle of the fiber (the relative change in density is negative).
The experimental time course data for the sarcomer length shown in Fig. 2 In the following analysis this measure is averaged over the defined central and peripheral  Fig. 2). It is worth mentioning that the amplitude of contraction or elongation of the fiber scales inversely with the fiber stiffness k: the softer the fiber, the stronger the mechanical deformation will be. Thus, a lower k value would simply explain the reported higher values for sarcomer contraction of about 30-40%. In our calculation we used k = 45nN/µm, a value reported by [48]. The experimentally measured equilibration time of the stress fiber upon stimulation is about 20 min (Fig. 2) which compares to about 3 min for the model results. These quantitave difference originate from two model simplifications.
First, we lump the focal adhesion associated processes into one equation thereby we shortcut the activation of ROCK and neglect prior activation steps of e.g. Rho-Gef or Rho-GTPase.
Considering these steps would cause an additional time delay. Secondly, the stimulation with calyculin happens instantaneously in the model omitting the time delay caused by the internalization of the drug. Refining the model and eliminating these simplifications will further decrease the differences in equilibration times.
To highlight the importance of the mechanical feedback we also include the expected results for a system neglecting this feedback shown as dashed lines in Fig. 11 and Fig. 12.
Here the homogeneous induction of the drug cause an almost uniform elevation of myosin activation within the cell. Slight differences between cell center and cell periphery would persist but rather marginally extend (Fig. 11). In fact stimulation leads also here to amplified distortions of the striation pattern. However, the changes in bandwidths are significantly smaller compared to the system incorporating the feedback (Fig. 12). Thus the closed biochemical and mechanical feedback loop is an essential feature required to describe the strong distortions of striation patterns upon homogenous drug induction. An intriguing aspect of our model is the way different stress fibers might cooperate inside a living cell. Conceptually it is easy to generalize our model to describe a system in which many stress fibers share the signaling input and many focal adhesions share the mechanical output. However, it remains a challenge to model also the dynamics of the actomyosin system if it is not completely condensed into stress fibers.

Model equations
Boundary conditions at x=0,L    FIG. 1: Cells adhere to the extracellular matrix by integrin-mediated contacts called focal adhesions. These contacts are the anchor points of stress fibers, which are actin filament bundles held together by the crosslinker molecule α-actinin and myosin II molecular motors. The myosins are assembled in myosin mini-filaments. Due to myosin motor activity stress fibers are under tension and exert forces to focal adhesions. This mechanical stimulus initiates biochemical signals (Rho-signal) that originate from focal adhesions and propagate into the cytoplasm, altering in turn myosin activity. Therefore the system of focal adhesions and stress fibers are connected by a biochemical and mechanical positive feedback loop (inset). The spatial part of our model is one-dimensional with one stress fibers extending between the two focal adhesions at x = 0 and x = L. Signaling pathway that controls myosin contractility depicted in its appropriate spatial context: mechanical cues are transduced to various biochemical signals at focal adhesions, however the precise mechanisms have not been resolved yet. One possible mechanism is that a Rho-GEF is activated by a mechanosensitive process at focal adhesions. Rho-GEF then promotes Rho-GTP loading and subsequent complexation with Rho-associated kinase (ROCK) which gets activated. Active ROCK is able to phosphorylate myosin light chain phosphatase (MLCP) at its myosinbinding subunit (MBS). MLCP and MLCP-P are freely diffusible in the cytoplasm and thus can reach the myosins in the stress fibers. Increased phosphorylation of MLCP to MLCP-P by ROCK effectively leads to increased phosphorylation of myosin light chain (MLC), thus increasing myosin contractility.
FIG. 4: Spatial dependence of the active myosin fraction n(x, t) at four different time points t ∈ {30s, 60s, 120s, 180s} as well as for the steady state, n ss (x): We solve the biochemical model implying the artificial initial conditions that all components are at zero activation level, but set the boundary force to 5nN . Because of this mechanical trigger at focal adhesions, MLC gets preferentially activated at the boundaries via the Rho-pathway which leads to a steady increase of the myosin activation level. Due to diffusible compounds in the Rho-pathway, the increased activation level is smoothed out towards the center of the cell. We model a stress fiber as a string of vec-elements such that the spring stiffness k n , the viscosity η n and the motor force F mn can vary spatially. E.g. the latter will vary spatially due to different myosin activity at the periphery compared to the center. The displacement of a certain site n is denoted by u n . The force F n onto this site n consists of elastic, viscous and motor contributions, compare Eq. (10). In the following we assume fix boundary conditions, namely u(0, t) ≡ 0 and u(L, t) ≡ 0 where L is the total length of the fiber.  14) which we include for comparison. For the assumed initial condition, ∂ x u(x, 0) = 0, the boundary force increases from its non-zero value at t = 0, given by Eq. (20) and quickly saturates at somehow larger values.  Fig. 4. For this simplifying case where the myosin activation level is not time dependent, Eq. (14) can be solved both numerically and analytically. Fig. 8 shows the analytical solution, Eq. (15), indicated by circles whereas the direct numerical solution of Eq. (14) is indicated by dotted lines. The solution is given at time points t ∈ {0.1s, 0.2s, 0.5s, 1.0s, 3.0s} as well as for the steady state u ss (x), assuming the initial condition: ∂ x u(x, 0) ≡ 0 and the boundary conditions: u(0, t) ≡ 0 and u(L, t) ≡ 0. . For the latter case, homogeneous induction of the drug cause an almost uniform elevation of myosin activity. Slight differences between center and periphery persist but rather marginally extend, whereas the closed feedback system result in an amplification of the spatial differences of myosin activation.