Cellular traction forces that are dependent on actin-myosin activity are necessary for numerous developmental and physiological processes [1]. Traction forces also contribute to undesired pathological states such as cancer metastasis [2–4]. Metastatic cancer cells generate larger traction forces, are more motile, and more invasive relative to their non-metastatic counterparts [3]. As traction force emerges as a promising cancer biomarker [5], there is a further need to understand force generation in response to chemical and mechanical cues. Here, with the help of computational, physics based models we aim to understand how much force a cell can really generate and what are the factors that limit this force. We are interested in the effects of four key components in traction force generation: number of active myosin motors, substrate stiffness, cell substrate adhesion dynamics, and mechanical reinforcement at adhesion complex.

The number of active motors can be controlled by biochemical signals that activate myosin light chain 2 (MLC2) phosphorylation. Number of active motors has been shown to alter force-velocity relationship [6–9]. While force generation is myosin dependent, force magnitude depends on other factors. Substrate stiffness has been shown to be an important parameter in cellular force generation. Stiffer substrates increase traction force [10–13]. Substrate stiffness also regulates the degree of cell-substrate adhesion and the size of the adhesion complex [14, 15]. Cells anchor to the substrate with adhesions ranging from 200 nm² to > 4 μm² [16, 17] adhesion size depending on the makeup of the adhesion plaque. Adhesions are protein plaques that have transmembrane components that interact with intercellular and extracellular components. Adhesion assembly is a hierarchical process [18]. Adhesion proteins are added as the plaque grows. The adhesome is variable between cell types but is highly complex and organized [18–20]. Regardless of the cell type, adhesion plaques physically connect the intracellular force-generating actin-myosin machinery to the cell’s substrate. This mechanical pathway allows cells to transfer the endogenous force to the substrate as traction force [21]. Of the proteins in the adhesion complex, integrins transmit endogenous force to the substrate [22–26]. Integrins span the membrane with an internal component where actin-myosin can attach via talin/vinculin and an external component that binds to substrate ligands [27–31]. Integrins are essential in force transmission as integrins have a force-dependent dissociation rate. As force increases, so does the probability that the integrin-substrate bond will rupture. Actin-myosin forces generated by the cell are transmitted as traction forces to the substrate via actin-myosin-integrin-substrate bonds.

We show stall force magnitude, the maximum force an actin-myosin contractile unit can generate and sustain without slipping, increases with actin-myosin activity and substrate stiffness. We also show that this max force is rarely achieved in contracting epithelial cells as the integrin bonds coupling actin-myosin to substrate break at far lower forces. Integrins with catch-slip bonds generate greater forces than slip bonds. Bond lifetimes for both slip and catch-slip integrins decrease with increasing actin-myosin activity and substrate stiffness. For these results, we focus on integrins attached so single actin-myosin filaments. Next, we probe how talin-mediated dynamics alter cellular force generation. Force dependent unravelling of talin can add additional actin-myosin filaments (recruitment) or cause talin branching to bind additional integrins (branching). Traction stress for catch-slip bonds with and without recruitment and branching exhibit biphasic or sigmoidal behavior that is dependent on the combination of myosin activity and substrate stiffness. Adhesion size, defined by the number of integrins actively contributing to traction stress, decreases with actin-myosin activity and substrate stiffness, however this is not a limiting factor in traction stress generation.

To understand how much force a cell can generate, we assume that traction stress, the force per unit area, is generated by actin-myosin filaments attached to integrin receptors [32–34]. We simulate a predefined adhesion plaque with a density of 100 integrins per 1 μm² [35, 36]. The adhesion plaque is a 10 × 10 grid representing 100 integrins precomplexed with talin and vinculin [37]. Talin is necessary to activate integrin [38–41] and does so by binding to integrin’s cytosolic domain [42]. Furthermore, talin’s force dependent unfolding is used to attach more filaments or integrins (Figure 1A). There are 100 actin-myosin filaments that are 3 μm [43] long, each containing 60 myosin motors [44–47]. Filaments bind to talin via vinculin [48–50].

The forces we are investigating originate from myosin binding and pulling along the actin filaments. We use a stochastic model of the myosin crossbridge cycle [51–53] modified with non-muscle myosin II transition rate parameters [54] to capture ATP hydrolysis, actin-myosin binding, conformational changes, and detachment (Figure 1B) [55, 56]. Simulations are performed using a timestep (Δt) of 1 ms. We model individual myosin motor head state changes stochastically, limiting individual motors to a single state change per timestep. Myosin binding sites on actin filaments are spaced by a fixed distance, bsd. Actin-myosin binding occurs when a myosin head is aligned to a binding site. Myosin undergoes thermal fluctuations that affects its ability to reach available binding sites. The distance is calculate based on the assumption that the potential energy of the myosin head is equal to the average unidirectional translational kinetic energy from the Boltzmann distribution, PE = ⟨KE⟩ where P E = 1 2 k m d 2 and ⟨ K E ⟩ = 1 2 k B T

The thermal fluctuation distance, d, that myosin can extend from its current position is ± d = k B T k m

where, k _m is the stiffness of the myosin motor stalk, defined below, k _B is the Boltzmann constant, and T is the temperature. Thus, a free ADP bound myosin head can bind to its respective actin filament if there is a binding site within this distance bsd ± d from myosin’s equilibrium position. The distance the myosin head has traveled away from it’s equilibrium unbound position to bind to the actin filament sets the initial strain, ϵ, in the myosin stalk, the coil-coiled region of the myosin tail that is not part of the thick filament [57]. The strain in the motor stalk, ϵ, is positive when strained in the direction that actin filament is actively sliding in along the myosin fibers. However, this strain in the myosin stalk/tail segment generates a force that opposes the sliding (Supplementary Figure S7).

Myosin state changes are calculated from the transition rates. The forward transition rates are k ₁₂, k ₂₃, k ₃₄, and k ₄₁ with k ₂₃ representing ATP hydrolysis and myosin head binding to actin. ATP hydrolysis and myosin binding to actin is a reversible process dependent on the ATP and ADP with k ₃₂ representing the reverse transition rate from state 3 to state 2. The crossbridge model consists of 4 states (state 1 through state 4) with respective transition rates from the current state, i, to the next state, j. Each motor state is independently updated during each time step. State transition rates follow the form k _ij .

Myosin dependent forces are regulated in part by the number of active myosin. We assume that only myosin with phosphorylated MLC2 can enter the crossbridge cycle and interact with the actin. Myosin that is not phosphorylated or that is dephosphorylated–phosphorylation being reversible–do not participate in the crossbridge cycle. Inactive myosins transition to state 15 from state 2 [58], modifying the 4-state cycle to a 5-state cycle [59]. Phosphorylation and dephosphorylation rates, k ₂₁₅ and k ₁₅₂ respectively, represent the transition rates from phosphorylated to dephosphorylated states. We define the ratio of phosphorylated MLC2 (ppMLC2) to non-phosphorylated MLC2 (MLC2) as k 152 k 215 = p p M L C 2 M L C 2 (1) where the phosphorylation ratio is the input parameter that modifies MLC2 phosphorylation.

Transition rates where the motor is bound to actin (states 1 and 3) are modified by strain in the motor’s stalk. The strain-dependent transition rate modifications are k s 12 = k 12 e 1 2 k m ϵ 2 k B T (2) and k s 34 = k 34 e k m ϵ δ k B T (3) where k _m is the stiffness in the current motor’s stalk, ϵ is the strain on the motor’s stalk, k _B is Boltzmann’s constant, T is the absolute temperature, δ is the characteristic bond length of the myosin bond [7], 1 2 k m ϵ 2 is the potential energy in the elastic motor head, and k _m ϵδ is the work done by the motor head. For now, the effect of strain on only these two transitions is considered as the rate of motor stepping is more responsive to forces during the first stage (3–4 transition), which is also the rate limiting step here, than the second stage (4–1 transition) [60].

Each motor’s state is independently updated. The probability, P _ij , that a motor will transition from state i to state j within a particular time step Δt is given by P i j t < Δ t = ∫ 0 Δ t k i j e − k i j t d t yielding P i j t < Δ t = 1 − e − k i j Δ t (4) for transition rates that are not modified by strain, and P i j t < Δ t = 1 − e − k s i j Δ t (5) for transition rates modified by strain. This relationship can also be arrived to by assuming the state transition reaction to be a binomial process with a fixed probability of success and a large number of trials.

States with only one transition rate (states 1, 15, and 4) transition to the next state (state 2, 2, and 1, respectively) when P _ij > rand, where rand is a uniformly distributed random number between 0 and 1. If P _ij ≤ rand, the motor stays in its current state. For states with two transition rates (states 2 and 3), states are updated based on joint probabilities. State 2 transitions as follows:

if P ₂₃ + P ₂₁₅ ≤ 1

and state 3 transitions as follows:

if P ₃₄ + P ₃₂ ≤ 1

Adhesion plaques form in a hierarchical process (Figure 1B [61], Zaidel-bar et al. [62]. The first step in forming an adhesion plaque requires integrin activation. Integrin activates by complexing with talin at a rate k _act . Activated integrins attach to the substrate ligand at a rate k _att . Once attached to the substrate, actin-myosin filaments bind to talin. These dynamics are temporally regulated. Integrin activation, integrin attachment, and actin-myosin attachment occur at independent time steps. The force generated in a single actin-myosin ensemble due to myosin activity is calculated by F a c t t = k m y N 1 + N 4 − ∑ i N k m ϵ i (6) where N ₁ and N ₄ are the number of motors in the ensemble in states 1 and 4 and y is the motor step size. The force in the actin-myosin filament is myosin-dependent. The first term represents the number of motors actively generating force by pulling on actin while the second term is the contribution of passive force from the strain of all actin-bound myosin in the ensemble. Negative strain in the motor stalk applies a force on the actin filament in the direction of active sliding between the actin filament and myosin bundles, while positive strain opposes active sliding (contraction). Motors that are not bound to actin do not contribute to the passive force.

Conversion of endogenous actin-myosin force to traction force requires a physical pathway. Integrins provide the physical link that determines if the endogenous force is converted to traction force. The force experienced at integrin-substrate bond is dependent on the displacement from pulling on integrin by the actin-myosin filament. We solve for the displacement F a c t − k s p r i n g x = γ d x d t (7) where x is the displacement, k _spring is the effective spring constant between substrate and integrin, and γ is the drag on the actin-myosin filament [63]. Integrin-substrate bond lifetimes have force-independent association rates and force dependent dissociation rates. Two integrin dissociation rates are reported: slip bonds and catch-slip bonds [64, 65]. Slip bond lifetimes continuously decrease with increasing force whereas catch-slip bonds have a force range that strengthens the bond and increases the bond lifetime [65, 66]. Catch-slip integrins experience an increase in bond lifetime up to an optimum applied force before reverting to a slip bond. Slip bond force-dependent dissociation rates are modeled by k o f f F i = k 0 e F i F b (8) where k ₀ is the unloaded off rate, F _i is the force at the integrin-substrate interface, and F _b is the characteristic bond rupture force. Integrins with catch-slip dynamics are modeled by [67]. k o f f F i = A e − F i ξ k B T + B e F i ξ k B T + C e − F i ξ k B T − 1 − 1 (9) where A, B, and C are constants, ξ is the unbinding length, and F _i is the load on the integrin. We use the dissociate rate, k _off , to calculate the dissociation probability P o f f = 1 − e − k o f f F i Δ t . (10)

Simulation for slip and catch-slip bond dynamics require talin-mediated integrin activation before integrin-substrate binding. The probability of integrin activating is P a c t i v a t e = 1 − e − k a c t Δ t (11) where k _act is the integrin activation rate. Activated integrins are then available to bind to the substrate. The probability for an integrin to bind to the substrate is P a t t a t c h t < Δ t = 1 − e − k a t t Δ t (12) where k _att it the rate at which activated integrin attach to the substrate. Activated and substrate-bound integrins are available for an actin-myosin contractile unit to attach. The probability that an actin-myosin filament attaches to integrin is P s f a t t t < Δ t = 1 − e − k s f a t t Δ t (13) where k _sfatt is the actin-myosin filament binding rate. We track integrin activation, integrin-substrate binding, actin-myosin filament attachment, and talin-mediated actin-myosin filament recruitment (more than one actin-myosin filament per integrin) or talin-integrin branching (more than one integrin per actin-myosin filament). We use talin-mediated dynamics to create a more realistic model. A complete physical pathway from actin-myosin to substrate enables transmission of endogenous forces to the cell’s immediate surrounding via integrin-substrate attachment. We calculate the force on individual integrin-substrate bonds by the actin-myosin filament activity as F i t = k s p r i n g x i t (14) where k _spring is the spring constant of the substrate and x _i is the displacement of the filament. We use this force in Eq. 9 to update the attachment status of integrins. In the case that the force does not rupture the integrin-substrate bond, the force is successfully transmitted to the substrate. Non-ruptured bonds allow myosin to continue to increase the displacement of the filament, increasing the forces experienced at the integrin-substrate bond. As force increases linearly with displacement there will be a max force at which the bond will break, ending its attachment period. Individual integrins experience force generated by actin-myosin. We use the individual forces at the integrins to determine the traction stress, or the collective force of all filaments actively generating and transmitting force. Traction stress is sensitive to bond lifetimes of individual integrins in the adhesion plaque. Traction stress is F T F = ∑ F i H i H i = 1  if physically linked H i = 0  if not physically linked . (15)

Similarly, adhesion size, which we define as only those integrins bound to the substrate actively transmitting actin-myosin force, is F A t = ∑ H i t H i = 1  if physically linked H i = 0  if not physically linked . (16)

For simulations where talin is permitted to recruit filaments or branch to integrins, we calculate the probability of talin unfolding and allow up to eight actin filaments to attach to a single integrin-talin complex or branching to two integrins. Talin unfolding probabilities are P r e c = 1 − e − k r e c Δ t (17) and P b r a n c h = 1 − e − k b r a n c h Δ t (18) with k _rec representing the talin unfolding rate for recruiting filaments and k _branch representing the talin unfolding rate for branching. We represent the adhesion plaque and the proteins required to establish a physical link from actin-myosin to substrate with matrices that are n × m where n is the number of integrins/filaments/attachment status, etc., and m is time in ms. Each matrix tracks motor state changes, integrin activation, integrin attachment, filament displacement, force, recruitment, and branching. For talin recruiting filaments, the displacement of the filaments is calculated by ∑ i = 1 R F a c t − k s p r i n g x = γ d x d t (19) where R is the number of recruited stress fibers. The force on the integrin is calculated using Eq. 14. Similarly the displacement of branched filaments attached to multiple integrins is calculated by F a c t − b k s p r i n g x = γ d x d t (20) where b is the number of branches. The force on the integrins is calculated using Eq. 14. Single filament average bond lifetimes are calculated by determining the time intact mechanical pathways are bound to the substrate τ ̄ = 1 n ∑ 1 n L t (21) where n is the number of intact pathways and L _t is the duration in ms of the binding event.

The stall force for the actin-myosin fibers is calculated by fitting the force generated over time with f = f s t a l l 1 − e − α t . (22)

The stall force reported is the average stall force over 100 actin-myosin filaments. When integrins bonds break in a force dependent manner, we obtain the maximum force generated by the actin-myosin filament during any attachment event. We report the average of the maximum force as F m a x a v g = 1 n ∑ i = 1 n F A M j (23) where F _AM is the vector n elements long that contains the maximum force generated in all binding events of filament j. Similarly, the average non-zero force generated is F n o n z e r o = 1 N ∑ F N Z (24) where N is the number of filaments and F _nz is the non-zero average of the filament. We report the maximum traction force per μm² by fitting F T F = F T F m a x 1 − e − β t . (25)

All model parameters and their values are given in Table 1.

The model described above explicitly simulates force and ATP dependent cycling of individual myosin motors as the enter and exit the crossbridge cycle and bind to and pull on the actin fibers. This allows us to exclude making any implicit assumptions regarding the force velocity dependence of myosin motors. The model also explicitly calculates the forces that these individual motors apply on actin fibers at each time point, how these forces are transmitted to the surface when the filament is bound to it via an elastic bond, and the dynamics of the individual bonds forming and breaking. In addition, the model simulates force dependent recruitment or branching of additional actin-myosin fibers at these bonds in order to simulate some of the mechanosensitive features of cell-substrate adhesion. We explore how nascent force develops in single filaments by varying integrin bond type, number of active motors, substrate stiffness, and talin-mediated actin-myosin filament recruiting or talin branching. Furthermore, we show how these factors influence traction stress transmitted by the actin-myosin fibers to the substrate.

The maximum force a single actin-myosin unit can generate is context-dependent and can be modulated by either one or more of the following conditions: integrin bond type, ppMLC2:MLC2, and substrate stiffness. Traction stress, the contribution of all force generating filaments in an adhesion complex, and force per filament have drastically different profiles even though traction stress is entirely dependent on single filament forces. Forces generated by actin-myosin filaments are modulated by two independent mechanisms: the number of motors and substrate stiffness. For example, simulations with integrin bonds that do not break show that if we hold ppMLC2:MLC2 constant, traction force increases with substrate stiffness. Likewise, if we keep substrate stiffness constant, stall force increases with ppMLC2:MLC2. Furthermore, force generation has a positive effect on myosin as myosin is mechanosensitive [74]. Myosin duty ratio shifts in a force-dependent manner [7]. The duty ratio increases with load which increases the number of attached motors as substrate stiffness increases in a ppMLC2:MLC2 independent mechanism. The increase in motors is not a response to stiffer substrates but a response to increased force. Stall force then is a response to both the number of activated motors and substrate stiffness. Our simulations with catch-slip bonds with and without talin-mediated actin-filament recruitment or integrin branching agree with experimental results showing that single integrin bonds withstand forces ranging from 10–30 pN [75, 76], up to > 54 pN [72]. Although catch-slip bonds can withstand forces greater than the 10–30 pN range reported, slip bonds rupture at forces < 30 pN (Figure 3A). As forces increase, bond lifetimes decrease. Interestingly, bond lifetimes for slip bond and catch-slip bonds settle at similar bond lifetimes while having different maximum force profiles. This is likely due to distribution of the bond lifetimes. At low ppMLC2:MLC2 and low substrate stiffness bond lifetimes for slip bonds skew farther right while at high ppMLC2:MLC2 and substrate stiffness catch-slip distributions skew farther right (Supplementary Figure S3). It is this extended positive skewness that increases the maximum force for catch-slip bonds. Bond lifetimes impact adhesion size. Increasing the force in individual actin-myosin filaments decreases the adhesion lifetime and consequently adhesion size. Adhesion size defined by only those integrins that are actively contributing to traction force is likely too narrow as it disregards scaffolding proteins that contribute to the composition of the adhesion. Nonetheless, our simulations show that depending on substrate stiffness, approximately 40–70% of the integrins are load bearing at any given time. These observation agree with experiments showing that approximately 60–75% of integrins in an adhesion are load bearing [71].

However, adhesion size seems to be a poor predictor for traction stress. Contrary to reports that use adhesions to predict traction force [17] there is a large range of traction stress at adhesions containing approximately 50 bound integrins (Figure 5A). Additionally, We see that traction stress is biphasic with adhesion size. There is an optimal adhesion size that maximizes traction force, past which, the traction force decreases. These results could explain how highly motile cells exhibit inverse correlation between adhesion size, migration speed, and increased invasiveness [3, 67, 77]. Increasing forces increases the likelihood that integrins will rupture their bond with the substrate and decreases the bond lifetime. Decreasing bond lifetimes affect adhesion size. We can make a generalization that large adhesion have integrins that experience longer bond lifetimes but transmit lower forces than small adhesions whose integrins experience shorter bond lifetimes but transmit greater forces. In this respect, our results agree with trends that relatively small adhesions transmit greater forces than large adhesions [78] and that adhesions < 1 μ m 2 exhibit widely variable traction stress [79, 80]. However, below a certain adhesion size, there might not be sufficient active attachments transmitting forces and traction force falls again. Individual integrins experience a range of forces that when evaluated as an aggregate reveal traction stress. Traction stress is not simply the product of motor activity, substrate stiffness, bond lifetime, or adhesion size. Experiments measuring traction stress report 3 trends: traction stress increases monotonically with substrate stiffness [81–88], increase then plateaus [89, 90], or decreases monotonically [91]. Our model predicts all 3 trends and under which conditions those trends occur. Within a specific substrate stiffness range spanning three orders of magnitude, regardless of ppMLC2:MLC2, we predict an increase in traction stress. Passing this stiffness threshold, traction stress plateaus for ppMLC2:MLC2 ≤ 0.2 or decreases then plateaus for ppMLC2:MLC2 ≥2. The traction stress biphasic behavior we observe is predicted in other models, either by fitting experimental data with a Gaussian curve [91] or by including actin retrograde flow, which occurs as actin-substrate anchoring proteins that act as molecular clutches disengage from the actin [92–94]. We note that the biphasic behavior of traction stress is context dependent. The biphasic behavior is not seen in simulations with low ppMLC2:MLC2. These simulations instead show a sigmoidal behavior. The existence of both biphasic and monophasic traction force response is a novel model prediction beyond what is reported in other models. While qualitative agreements between the model predictions and various experimental observations can be made, further testing of specific model outcomes is required to validate the model assumptions. For example, manipulating myosin activation rates by using Rho/ROCK blockers and agonists, integrated with single molecule force spectroscopy in cells attached to or embedded within substrates with different stiffness would further help test a number of specific model predictions made above.

Overall, the model presented here is an extension of, or borrows from, a number of prior actin-myosin force generation and filament sliding models [51, 57, 95]. It also borrows from other traction force, actin-integrin-substrate adhesion models [13, 36, 96]. Specific features new to this model are—1) The model describes five stages to allow for myosin activation via phosphorylation, as compared to the one [51]. This allows the future possibility of integrating it with systems biology models focused on upstream signals that can promote or inhibit myosin phosphorylation. Additionally, it incorporates mechanical regulation of cross-bridge cycling rates as a function of individual motor strains, fluctuations in myosin distributions along the actin-myosin fibers, and the ability to track activity at the level of individual motors. Beyond the detailed incorporation of actin-myosin cross-bridge cycling, the model also incorporates differences in cell-substrate adhesion bond kinetics, substrate mechanical stiffness, mechanochemical regulation of focal adhesion formation via Talin binding and mechanosensitive recruiting or branching of actin stress-fibers, along the lines of the model described by [63]. These allow for easy extension of the model to probe biochemical signaling downstream of mechanosensitive focal adhesions. While we limit the use of the model to focus on factors that govern instantaneous traction force generation in cells, there are many other potential applications of this model, especially when coupled with biochemical signaling models. A limitation of this model is that it tracks the generation of traction forces only on the order of minutes. Thus, cell scale events such as cytoskeletal remodeling ([97], actin fiber depolymerization and reorganization [98], adhesion remodeling occurring over longer time-scales [99, 100], as well as changes to protein expression levels within the cell cannot be accounted for. These phenomena could transiently alter traction stress which might differ from the short timescale steady state predictions made here. The bridging of phenomena across times scales and observing transient changes to traction stress generation is a challenge that we look to addressing in the future iterations and applications of this model.

Different cell types generate different levels of traction stress. The forces themselves are not exclusive to a cell type. Rather, forces are dependent on actin-myosin organization, myosin contractility, cell-substrate adhesion and coupling between the cell cytoskeleton and the cell-substrate adhesions. These factors are in-turn governed by intra-cellular and extra-cellular chemical and mechanical signals such as hypoxia, adrenergic signaling, growth factor signaling, substrate viscoelasticity and so on. [82, 101, 102]. In this work, we focused on the immediate regulators of cell traction force generation—myosin motor activity (primarily non-muscle myosin 2 interactions with actin filaments), actin-myosin filament coupling to cell-substrate adhesions, and substrate elastic response. We find that stall forces are not the limiting factor for the traction force that cells can generate. Instead, as substrate stiffness increases the number of active motors in the force generating states increase, consequently increasing the stall force for individual actin-myosin filaments. Investigating the actin-integrin-substrate adhesion dynamics, we find that traction stress is limited by the integrin-substrate bond lifetime, which decreases with increasing substrate stiffness and myosin activity. Ultimately, we show that traction stress is proportional to the impulse transmitted by individual actin-myosin filaments to the substrate via focal adhesions. Our model integrates the nanoscale, microsecond activity of individual myosin motors interacting with actin filaments with nanoscale cell substrate adhesion mechanics to estimate traction forces generated by cells over microscale adhesion regions over a timescale of minutes. It also allows for ways to integrate biochemical signals that alter myosin phosphorylation states, actin organization, focal adhesion protein concentrations and activity, as well as mechanochemical feedback between these elements to provide a detailed cell traction force generation model. The results of this model can further be integrated with systems biology models focusing on intra-cellular chemical signaling to include the effects of mechanical force generation on certain mechanically activated reaction pathways such as talin-vinculin binding and focal adhesion kinase (FAK) activation and signaling [103].