Disordered actomyosin networks are sufficient to produce cooperative and telescopic contractility

While the molecular interactions between individual myosin motors and F-actin are well established, the relationship between F-actin organization and actomyosin forces remains poorly understood. Here we explore the accumulation of myosin-induced stresses within a two-dimensional biomimetic model of the disordered actomyosin cytoskeleton, where myosin activity is controlled spatiotemporally using light. By controlling the geometry and the duration of myosin activation, we show that contraction of disordered actin networks is highly cooperative, telescopic with the activation size, and capable of generating non-uniform patterns of mechanical stress. We quantitatively reproduce these collective biomimetic properties using an isotropic active gel model of the actomyosin cytoskeleton, and explore the physical origins of telescopic contractility in disordered networks using agent-based simulations.


Supplementary Tables
Supplementary Table 1

Supplementary Note 1: Linear Active Gel Model
We describe a continuum model for the dynamics of disordered actomyosin based on active gel theories [15]. We assume that the thin sheet of actomyosin gel is in mechanical equilibrium with the substrate and the internal mechanics of the gel is described by a stress tensor, , where i, j denote in-plane spatial coordinates. In the absence of adhesion, we assume that the substrate provides a viscous drag to the gel characterized by a friction tensor . In-plane force balance in the thin film limit gives us, , where h is the gel thickness. Here we suggest a model for the mechanics of the gel by defining the constitutive relations for material stress.
The total stress is decomposed as the sum of elastic, dissipative and active stresses, . While elasticity arises primarily from the compliance of F-actin filaments, dissipative stresses can arise from processes such as disentanglement of filaments and myosin binding/unbinding. Using the plane stress approximation (in the thin-film limit), the constitutive relations for elastic and dissipative stresses are given by, where is the symmetrized strain tensor, is the strain-rate tensor, E is the Young's modulus, is the poisson ratio, and are the bulk and shear viscosities. For simplicity we assume that myosin induced active stresses are isotropic, with the functional form, where is the magnitude of contractile stress, is a positive constant (Hill coefficient), is the myosin density (assumed to be uniform), is the density beyond which myosin induced stresses saturate, and is the timescale for the accumulation of active stresses. To model the experimental geometry, we assume a spatially constant active stress profile for (and zero outside), where is the radius of the activation zone. We now exploit radial symmetry of the flow profile as observed experimentally to assume that all the quantities are azimuthally invariant. In polar coordinates, force balance in the radial direction is given by the stress equilibrium condition (Eq. S17). In terms of the radial displacement field ( ), the force balance equation simplifies to, , where, is the viscoelastic timescale and The solution to (S6) can be obtained analytically using the initial condition of zero displacement (the solution is cumbersome and hence not presented here). Using the analytical solution for radial displacement we can then evaluate mechanical quantities such as radial strain ( ), orthoradial strain ( ), velocity ( ), strain-rates ( and ) and the corresponding stresses.

Supplementary Note 2: Active Gel Model Myosin Dynamics
In our continuum model we treat actomyosin as a single unit, and assume a uniform active stress profile. The active stresses arise from myosin that are highly dynamic. The concentration dynamics of myosin minifilaments (average length l) bound to actin can be modeled using the following equation [16]: where D is a diffusivity consant, is the unbinding rate of myosin, is the binding rate and is the total concentration of myosin. The equation tells us that density fluctuations in myosin relax over a characteristic timescale, 1/ / , beyond which myosin is primarily advected by actin flow.

Supplementary Note 3: Agent-Based Simulation
Displacement of each cylindrical segment is governed by the Langevin equation with inertia neglected: where i r is the position of endpoint of cylindrical segments, ζ i is an effective drag coefficient, t is time, and i F is a net deterministic force including extensional, bending, and repulsive forces.
T i F is a stochastic force determined by the fluctuation-dissipation theorem [6]: where δ ij is the Kronecker delta, δ is a unit second-order tensor, and Δt is time step. For the cylindrical geometry of segments, we used the approximate form of ζ i [7]: where c,i r and 0,i r are the diameter and length of a segment, respectively. Position of each segment is updated using the Euler integration scheme: The extension and bending of actin and motor are governed by harmonic potentials with stiffness s  and b  , respectively: where r is a distance, θ is a bending angle, and the subscript 0 denotes the equilibrium value. Repulsive forces that simulate volume-exclusion effects are calculated via the following harmonic potential [9]: where r  is strength of repulsive force, and 12 r is the minimum distance between a pair of actin segments.
Each of the motor arms attached to the backbone represents 8 myosin heads (N h = 8). Thus, each motor with 8 arms behaves as a thick filament with 64 myosin heads. A free arm can bind to a binding site on F-actin at a rate, 40·N h ·s -1 . Motor arms walk on and unbind from F-actin following the walking ( ) and unbinding rates ( ) which vary depending on forces applied to the arms (Supplementary Fig. 8). At each walking event, arms slide from a current binding site to a next one located toward the barbed end by ~7 nm. After reaching the barbed end, motors slide off from F-actin via a next walking event. Note that it is assumed that myosin heads behave as a catch bond [10,11], leading to lower and with higher applied forces. We assume that deactivated motors located outside the activation zone unbind at a much higher rate (

Supplementary Note 4: Agent-Based Model Abstraction
Blebbistatin binds the myosin-ADP-Pi complex, and inhibits phosphate release. Thus, it keeps myosin in the actin-detached state, preventing rigid cross-linking [13,14] effectively reducing processivity. To this effect, in our simulations, we increase myosin unbinding in regions where blebbistatin is active and decrease myosin unbinding in regions where blebbistatin has been inactivated. In doing so, the F-actin outside the activation region is loosely crosslinked, and thus as contraction proceeds in the activation region, the F-actin separates between these two regions.

Supplementary Note 5: Nonlinear Active Gel Model
In the nonlinear active gel model we assume that the elastic modulus and the viscosity are functions of actomyosin density, , such that where and are positive constants. As a consequence, the network stiffens and becomes more viscous with increasing density. The constitutive laws for elastic and viscous stresses are assumed to be linear in strain and strain rates, respectively, as defined in Eqs. (S2) and (S3).
Since changes in density are slaved to the network strain, we have the continuity equation: . (S16) Using this simple form for strain-dependent viscoelastic response, we numerically solve the resultant force balance equations (Mathematica, NDSolve), Eq. (S1), assuming spatially uniform contractile stress within the activation domain. By varying the parameters and in the model, we can control the strength of nonlinear coupling and density dependent variations in mechanical properties. We find that the telescopic behavior is retained for non-zero values of and ( Supplementary Fig. 9A), but the slope of the relationship between velocity and activation is reduced for networks stiffening ( 0 and thickening ( 0 with increasing density. This attenuation in telescopic contractility results from concomitant reductions in sptially averaged maximum strain and strain-rates ( Supplementary Fig. 12 B,C).

Velocity Calculations
The velocity field of the F-actin is calculated using Particle Image Velocimetry (PIV) (mPIV; www.oceanwave.jp/softwares/mpiv/). mPIV creates a regular grid across the image, where the velocity is calculated for each grid location between a pair of images. The strain is taken as the divergence of the F-actin displacement field, and the strain rate is the divergence of the velocity field, for the experiment, model and simulation. For each metric, the average is taken over the activation zone. This leads to strains larger than 2, as in addition to the sum of radial and orthoradial components, there is additional F-actin which flows into the activation region over the time at which the strain and strain rate are averaged. The velocity is averaged only within a narrow vicinity at the boundary, not throughout the activation area.

Mechanical Stress Calculation
In mechanical equilibrium, Eq. (S1) provides the condition for local force-balance in the gel. Now assuming radial symmetry, mechanical equilibrium in polar coordinates satisfies the where and are the two normal stresses in the radial and the orthoradial directions. The general solution can be expressed in terms of an Airy stress function [4], . (S19) Since is a biharmonic function, we can use the ansatz [5], .
Using the condition that is finite at r=0, we have b=0 and c=0. To determine the constant a, we use stress boundary condition 0, where R defines the total system size ( ≫ ). The profile for radial stress can then be obtained from the experimentally measured flow profile, up to an undetermined multiplicative factor, , . (S20) Since the flow velocity points radially inward ( 0) and localizes at the boundary of the activation region, the internal stress is positive and largest at the center of the activation region where it accumulates over time (Fig. 2E).

Principal Component Analysis
To analyze the local shape and magnitudes of strain and strain-rates we perform a principal component analysis of the symmetric strain tensor [4], defined as: where u is the displacement field. The symmetric strain-rate tensor is given by the time derivative of the strain tensor, or equivalently as the spatial gradients of the velocity field. The principal strains, s 1 and s 2 , are defined by the eigenvalues of the strain tensor. Their magnitudes define the major and the minor axes of the principal strain ellipse, which characterizes the shape of local strains. For example, isotropic strains correspond to , whereas uniaxial deformations imply 0 or 0. The orientations of the ellipse axes correspond to the eigenvectors of the strain tensor, such that the major (minor) axis points along the direction of maximum (minimum) normal strain. Total normal strain is given by , which equals the divergence of the displacement field. The local maximum value of the shear strain, , is related to the difference of the principal strains: | |/2 [4].

Hill Model Regression
The regression to the Hill Model ( Fig 3C) Three parameters, were fit, V max , K 0.5 and n H are as follows: 2.48, 0.56 m -2 , and 10.95 respectively. The Hill coefficient represents the cooperativity of myosin behavior towards generating strain. More specifically, assuming that the density of myosin thick filaments is proportional to the total density of actomyosin bonds, the hill coefficient represents the sensitivity of the strain to the number of actomyosin bonds. V max and K 0.5 represent the maximum strains accomplished during contraction, and the quantity of myosin it would take to reach half of that strain.

Orientation Quantification
Analysis of the orientation of F-actin filaments was performed using the ImageJ (NIH) plugin OrientationJ (EPFL; http://bigwww.epfl.ch/demo/orientation/). The color is orientation, the saturation is coherence, and the brightness is the original image brightness. The global coherency of F-actin is measured by and represents the ratio between the difference and the sum of the tensor eigenvalues, and is bounded between 0 (isotropic) and 1 (aligned) [1].

Order Parameter Calculation
The nematic order parameter, q was calculated as performed previously [2] [3] using Custom written Matlab (Mathworks, Natick, MA) routines. Briefly, the image was broken down into a series of small windows (3.5 x 3.5 µm) to determine the local orientation of filaments. For each window, a symmetric 2D Gaussian filter was applied, and then the 2D FFT was determined. The orientation of the FFT transform was determined by calculating the least second moment of rotation. The orientation in real space is orthogonal to the orientation calculated in frequency space. To calculate the order parameter, q, a given orientation vector was compared to its nearest neighbors using the following equation: The order parameter thus calculates the degree of F-actin alignment, with 0 being unaligned, and 1 being highly aligned.