article: Actomyosin-based cortical flow is a fundamental engine for cellular morphogenesis. Cortical flows are generated by cross-linked networks of actin filaments and myosin motors, in which active stress produced by motor activity is opposed by passive resistance to network deformation. Continuous flow requires local remodeling through crosslink unbinding and and/or filament disassembly. But how local remodeling tunes stress production and dissipation, and how this in turn shapes long range flow, remains poorly understood. Here, we study a computational model for a cross-linked network with active motors based on minimal requirements for production and dissipation of contractile stress: Asymmetric filament compliance, spatial heterogeneity of motor activity, reversible cross-links and filament turnover. We characterize how the production and dissipation of network stress depend, individually, on cross-link dynamics and filament turnover, and how these dependencies combine to determine overall rates of cortical flow. Our analysis predicts that filament turnover is required to maintain active stress against external resistance and steady state flow in response to external stress. Steady state stress increases with filament lifetime up to a characteristic time τm, then decreases with lifetime above τm. Effective viscosity increases with filament lifetime up to a characteristic time τc, and then becomes independent of filament lifetime and sharply dependent on crosslink dynamics. These individual dependencies of active stress and effective viscosity define multiple regimes of steady state flow. In particular our model predicts that when filament lifetimes are shorter than both τc and τm, the dependencies of effective viscosity and steady state stress on filament turnover cancel one another, such that flow speed is insensitive to filament turnover, and shows a simple dependence on motor activity and crosslink dynamics. These results provide a framework for understanding how animal cells tune cortical flow through local control of network remodeling.
  Cortical flow is a fundamental and ubiquitous form of cellular deformation that underlies cell polarization, cell division, cell crawling and multicellular tissue morphogenesis [1–6]. Cortical flows originate within a thin layer of cross-linked actin filaments and myosin motors, called the actomyosin cortex, that lies just beneath the plasma membrane [7]. Local forces produced by bipolar myosin filaments are integrated within cross-linked networks to build macroscopic contractile stress [8–10]. At the same time, cross-linked networks resist deformation and this resistance must be dissipated by network remodeling to allow macroscopic deformation and flow. How force production and dissipation depend on motor activity and network remodeling remains poorly understood. One successful approach to modeling cortical flow has relied on coarse-grained phenomenological descriptions of actomyosin networks as active fluids, whose motions are driven by gradients of active contractile stress and opposed by an effectively viscous resistance [11]. In these models, spatial variation in active stress is typically assumed to reflect spatial variation in motor activity and force transmission [12], while effective viscosity is assumed to reflect the internal dissipation of elastic resistance due to local remodeling of filaments and/or cross-links [7,13]. Models combining an active fluid description with simple kinetics for network assembly and disassembly, can successfully reproduce the spatiotemporal dynamics of cortical flow observed during polarization [11], cell division [14,15], cell motility [16,17] and tissue morphogenesis [18]. However, it remains a challenge to connect this coarse-grained description of cortical flow to the microscopic origins of force generation and dissipation within cross-linked actomyosin networks. Studies in living cells reveal fluid-like stress relaxation on timescales of 10-100s [1,2, 11,19–21], which is thought to arise through a combination of cross link unbinding and actin filament turnover [7,13,22]. Theoretical [23,24] and computational [25–27] studies reveal that cross-link unbinding can endow actin networks with complex time-dependent viscoelasticity. However, while cross-link unbinding is sufficient for viscous relaxation (creep) on very long timescales in vitro, it is unlikely to account for the rapid cortical deformation and flow observed in living cells [26,28–31]. Experimental studies in living cells reveal rapid turnover of cortical actin filaments on timescales comparable to stress relaxation (10-100s) [32–35]. Perturbing turnover can lead to changes in cortical mechanics and in the rates and patterns of cortical flow [33,36]. However, the specific contributions of actin turnover to stress relaxation and how these depend on network architecture remain unclear. Recent work has also begun to reveal mechanisms for active stress generation in disordered actomyosin networks. Theoretical studies suggest that spatial heterogeneity in motor activity along individual filaments, and asymmetrical filament compliance (stiffer in extension than in compression), are sufficient for macroscopic contraction [37,38], although other routes to contractility may also exist [38]. Local interactions among actin filaments and myosin motors are sufficient to drive macroscopic contraction of disordered networks in vitro [39], and the kinematics of contraction observed in these studies support a mechanism based on asymmetrical filament compliance and filament buckling. However, in these studies, the filaments were preassembled and network contraction was transient, because of irreversible network collapse [40], or buildup of elastic resistance [41], or because network rearrangements (polarity sorting) dissipate the potential to generate contractile force [42–45]. This suggests that network turnover may play essential role (s) in allowing sustained production of contractile force. Recent theoretical and modeling studies have begun to explore how this might work [46–48], and to explore dynamic behaviors that can emerge when contractile material undergoes turnover [15,49]. However, it remains a challenge to understand how force production and dissipation depend individually on the local interplay of network architecture, motor activity and filament turnover, and how these dependencies combine to mediate tunable control of long range cortical flow. Here, we construct and analyze a simple computational model that bridges between the microscopic description of cross-linked actomyosin networks and the coarse grained description of an active fluid. We represent actin filaments as simple springs with asymmetric compliance; we represent dynamic binding/unbinding of elastic cross-links as molecular friction [50–52] at filament crossover points; we represent motor activity as force coupling on a subset of filament cross-over points with a simple linear force/velocity relationship [53]. Finally, we model filament turnover by allowing entire filaments to appear with a fixed probability per unit area and disappear with fixed probabilities per unit time. We use this model to characterize: first the passive response of a cross-linked network to externally applied stress, then the buildup and maintenance of active stress against an external resistance, and finally the steady state flows produced by an asymmetric distribution of active motors in which active stress and passive resistance are dynamically balanced across the network. Our results reveal how network remodeling can tune cortical flow through simultaneous effects on active force generation and passive resistance to network deformation.
  We assume (Fig 1A) that local deformation of filament i gives rise to equal and opposite elastic restoring forces on its endpoints: F p, i elas = μ γ i u i ^, F b, i elas = - F p, i elas (1) where γi = (|bi − pi| − L) /L is the strain on filament i, and μ is a normalized spring constant. To model asymmetric filament compliance, we use a piecewise linear approximation to the non-linear entropic force extension curve for semi-flexible polymers with lengths less than the persistence length [56,57]. We set μ = μe if the strain is positive (extension), and μ = μc if the strain is negative (compression). Previous models have represented cross-linkers as elastic connections between pairs of points on neighboring filaments that bind and unbind with either fixed or force-dependent probabilities [25,58]. Here, we introduce a coarse-grained representation of crosslink dynamics by introducing an effective drag force that couples every pair of overlapping filaments, and which represents a molecular friction arising from the time-averaged contributions of many individual transient crosslinks (Fig 1B). This is a reasonable approximation if all pairs of overlapping filaments have equal access to a non-limiting pool of cross links, and if the rate at which filaments move past one another is slow relative to the unbinding rate of individual crosslinks [59]. This coarse-grained approach has been used to model frictional forces arising from ionic cross-linking of actin filaments in vitro [60,61], and simple force-velocity relationships for systems of cytoskeletal filaments and cross-linking motors [53,62–64]. To implement coupling through effective drag, for any pair of overlapping filaments i and j, we write the drag force on filament i as: F i, j ξ = - ξ (v i - v j) (2) where ξ is the drag coefficient and vi, vj are the average centroid velocities of filaments i and j. We apportion this drag force to the two endpoints pi and bi of filament i as follows: If xi, j is the position of the filament overlap, then we define λi, j = |xi, j − pi|/|bi − pi| to be the fractional position of the overlap point along filament i, and we assign (1 - λ i, j) F i, j ξ to endpoint pi and λ i, j F i, j ξ to endpoint bi. The total crosslink coupling forces on endpoints pi and bi, due to overlaps along filament i, can then be written: F p, i xl = ∑ j (1 - λ i, j) F i, j ξ F b, i xl = ∑ j λ i, j F i, j ξ (3) where the sums are taken over all filaments j that overlap with filament i. This model assumes a linear relation between the drag force and the velocity difference between attached filaments. Although non-linearities can arise through force dependent detachment kinetics and/or non-linear force extension of cross-links, we assume here that these non-linear effects are of second or higher order. We add motor activity at the point of overlap between two filaments i and j as follows: For each filament in the pair, we impose an additional force of magnitude υ, directed towards its pointed (-) end (Fig 1C): F i υ = - υ u i ^ (4) We impose an equal and opposite force on its overlapping partner. We distribute these forces to filament endpoints as described above for crosslink coupling forces. Thus, the total force on endpoints i and i+1 due to motor activity on overlap points between filaments i and j can be written as: F p, i motor = ∑ j (1 - λ i, j) (F i υ - F j υ) q i, j = υ ∑ j (1 - λ i, j) (u j ^ - u i ^) q i, j F b, i motor = υ ∑ j (λ i, j) (u j ^ - u i ^) q i, j (5) where j enumerates over all filaments j that overlap with filament i, and qi, j equals 0 or 1 depending on whether there is an “active” motor at this location. To model dispersion of motor activity, we set qi, j = 1 on a randomly selected subset of filament overlaps, such that q ¯ = ϕ, where q ¯ indicates the mean of q (Fig 1C). To write the full equation of motion for a network of actively and passively coupled elastic filaments, we assume the low Reynold’s number limit in which inertial forces can be neglected, and we equate the sum of all forces acting on all filament endpoints to zero to obtain: 0 = - ζ v i p L / 2 - F p, i xl + F p, i elas + F p, i motor 0 = - ζ v i b L / 2 - F b, i xl + F b, i elas + F b, i motor (6) where the first terms in each equation represent the hydrodynamic drag on the half-filaments adjoining endpoints pi or bi with respect to motion at velocities v i p or v i b against the surrounding fluid, and ζ is the drag coefficient. We used a mikado model approach [65] to initialize a minimal network of overlapping unstressed linear filaments in a rectangular 2D domain. We generate individual filaments by laying down straight lines, of length L, with random position and orientation. We define the density using the average distance between cross-links along a filament, lc. A simple geometrical argument can then be used to derive the number of filaments filling a domain as a function of L and lc [66]. Here, we use the approximation that the number of filaments needed to tile a rectangular domain of size Dx × Dy is 2DxDy/Llc, and that the length density is therefore simply, 2/lc. Although we do not model thermal forces explicitly, the contribution of thermal fluctuations to filament elasticity are embedded in our coarse-grained representation of asymmetrical filament compliance. In principle, thermally-driven transverse fluctuations of filament segments between crosslink points could influence crosslink binding kinetics. However, for the network mesh sizes considered here, lc <= 0. 5μm, the root mean square amplitude of these fluctuations is predicted to be < 5nm (see e. g [56]), suggesting that these effects will be minor. Hence, we have chosen to ignore them here. In living cells, actin filament assembly is governed by multiple factors that control filament nucleation, branching and elongation. Likewise filament disassembly is governed by multiple factors that promote filament severing and monomer dissociation at filament ends. Here, we implement a very simple model for filament turnover in which entire filaments appear with a fixed rate per unit area, kapp and disappear at a rate kdissρ, where ρ is the filament density (Fig 1D). With this assumption, in the absence of network deformation, the density of filaments will equilibrate to a steady state density, kapp/kdiss, with time constant τr = 1/kdiss. In deforming networks, filament density will also decrease under extensional strain and increase under compressional strain. Thus filament density will be set by a dynamic interplay of deformation and density equilibration via turnover (see below and (S1 Appendix A. 3) ). To implement this model, at fixed time intervals τs < 0. 01 ⋅ τr (i. e. 1% of the equilibration time), we selected a fraction, τs/τr, of existing filaments (i. e. less than 1% of the total filaments) for degradation. We then generated a fixed number of new unstrained filaments kappτsDxDy at random positions and orientations within the original domain. We refer to kdiss = 1/τr as the turnover rate, and to τr as the turnover time. Further details regarding our simulation approach and references to our code can be found in the Supplementary Information (S1 Appendix A. 1). Briefly, eqs 1–6 define a coupled system of ordinary differential equations that can be written in the form: A · x ˙ = f (x) (7) where x is a vector of filament endpoint positions, x ˙ the endpoint velocities, A is a matrix with constant coefficients that represent crosslink coupling forces between overlapping filaments, and f (x) represents the active (motor) and elastic forces on filament endpoints. We smoothed all filament interactions, force fields, and constraints linearly over small regions such that the equations contained no sharp discontinuities. We used a fourth-order Runga-Kutta method to numerically integrate this system of equations to find the time evolution of the positions of all filament endpoints. We generate a network of filaments with random positions and orientations as described above within a domain of size Dx by Dy. For all simulations, we imposed periodic boundaries in the y-dimension. To impose an extensional force per unit length (2D stress) on the network, we constrained all filament endpoints within a fixed distance 0. 05 ⋅ Dx from the left edge of the domain to be non-moving, then we imposed a rightwards force on all endpoints within a distance 0. 05 ⋅ Dx from the right edge of the domain, such that the force per unit length of boundary equals the desired stress value. To simulate free contraction, we removed all constraints at domain boundaries; to assess buildup and maintenance of contractile stress under isometric conditions, we used periodic boundary conditions in both x and y dimensions. In our 2D model, we measure stress as a force per unit length. We measured the internal network stress at each axial position by summing the axial (x) component of the tensions on all filaments intersecting that position, and dividing by network height Dy. We quantified two different forms of strain: the average filament strain, which measures the deformations of individual filaments, and the cumulative network strain, a normalized measure of network deformation defined as the change in axial length of a patch of network divided by its original length. These two measures can differ because filaments can slide relative to one another during deformation and because strained filaments are replaced by unstrained filaments during network turnover. We measured the strain on individual filaments as defined above from γi = (|bi − pi| − L) /L. Then we averaged this measurement over all filaments in a network to obtain an average filament strain. To measure the average network strain rate, we first measured the mean velocity v (X) at position X (relative to the network boundary at x = 0) to be the average velocities of all filaments intersecting that position. In the cases where we measure network strain or strain rates, we observed an approximately linear dependence of v (X) on X; hence the strain rate is approximately uniform across the network (Fig 2B and S6 Fig). Accordingly, for each filament, we took 1 X d X d t to be an estimate of the strain rate on the network between x = 0 and x = X. We averaged this estimate over all filaments in a domain to get an average strain rate. Finally, to estimate the cumulative network strain at a given time T in the simulation, we integrated the strain rate with respect to time for t = 0 to T. We assigned biological plausible reference values for all parameters (See Table 1). For individual analyses, we sampled the ranges of parameter values around these reference values shown in S1 Table.
  Thus far, we have considered independently how network remodeling controls the passive response to an external stress, or the steady state stress produced by active contraction against an external resistance. We now consider how these two forms of dependence can combine to shape steady state flow produced by spatial gradients of motor activity. To this end, we model a simple scenario in which motor activity is continuously patterned such that the right half network has uniformly high levels of motor activity (controlled by υ, with ϕ = 0. 5), while the left half network has none (ϕ = 0). For simplicity, we imposed periodic boundary conditions at left and right boundaries. Under these conditions, with filament turnover, we expect the right half network to contract continuously against a passive resistance from the left half network. Given the highly asymmetric filament compliance, the internal resistance of the right half network to active compression should be negligible compared to the external resistance of the left half network to extension. Thus the steady state flow should be well-described by: γ ˙ = σ s s η (11) where σss is the active stress generated by the right half-network (less the internal resistance to filament compression), η is the effective viscosity of the left half network and strain rate γ ˙ is measured in the left half-network. Note that strain rate can be related to the steady state flow velocity v at the boundary between right and left halves through v = γ ˙ D x. Therefore, we can understand the dependence of flow speed on filament turnover and other parameters using the approximate relationships summarized by eqs 9 and 10 for η and σss. As shown in Fig 8, there are two qualitatively distinct possibilities for the dependence of strain rate on τr, depending on the relative magnitudes of τm and τc. In both cases, for fast enough turnover (τr < min (τm, τc) ), we expect weak dependence of strain rate on τr (γ ˙ ∼ τ r 1 / 4). For all parameter values that we sampled in this study (which were chosen to lie in a physiological range), τm > τc. Therefore we predict the dependence of steady state strain rate on τr shown in Fig 8A. To test this prediction, we simulated the simple scenario described above for different values of τr, with all other parameter values initially fixed at their reference values (Fig 9A). As expected, for all values of τr, the asymmetric pattern of active contraction gave rise to steady state flow, characterized by continuous contraction of the right half-network and expansion of the left half-network, with peak velocity at the boundaries between right (contracting) and left (expanding) domains (Fig 9B). At long times, the average strain on individual filaments reached a plateau (see S6 Fig), but the cumulative network strain increased linearly with time (Fig 9C), indicating steady state flow with a constant strain rate and peak velocity. Plotting steady state strain rate vs filament lifetime τr confirmed the predicted dependence: Steady state strain rates approached zero with increasing τr; however, for decreasing τr, steady state strain rates increased steadily, before reaching an approximate plateau on which strain rate varied by less than 15% over more than two decades of variation in τr (Fig 9D). We repeated these simulations for a wider range of parameter values, and saw the same qualitative dependence of γ ˙ on τr in all cases. Using eq 9 with τr < τc and eq 10 with τr < τm, and the theoretical or empirical scaling relationships found above for ηc, τc, σm and τm, we predict a simple scaling relationship for γ ˙ (for small τr, see Fig 9D): γ˙=υξL (τr) 1/4 (12) Indeed, when we plot the steady state measurements of γ ˙, normalized by υ/ξL, for all parameter values, the data collapse onto a single curve for small τr. Thus. our simulations identify a flow regime, characterized by sufficiently fast filament turnover, in which the steady state flow speed is buffered against variation in turnover, and has a relatively simple dependence on other network parameters.
  Cortical flows arise through a dynamic interplay of force production and dissipation within cross-linked actomyosin networks. Here we combined computer simulations with simple theoretical analysis to explore how this interplay depends on motor activity, crosslink dynamics, network architecture and filament turnover. Our results reveal two essential requirements for filament turnover during cortical flow: (a) to allow the continuous relaxation of elastic resistance without catastrophic loss of network connectivity and (b) to prevent the dissipation of active stress through local network rearrangements. We find that biphasic dependencies of active stress and passive relaxation on filament lifetime define multiple modes of steady state flow with distinct dependencies on network parameters and filament turnover. We identify two distinct modes of passive response to uniaxial stress: a low turnover mode in which filaments strain to an elastic limit before turning over, and effective viscosity depends on crosslink density and effective crosslink friction, and a high turnover mode in which filaments turn over before reaching an elastic limit and effective viscosity is proportional to elastic resistance and approximately proportional to filament lifetime. We note that the weakly sub-linear dependence of effective viscosity on filament lifetime that we observe in the high turnover regime may simply reflect a failure to capture very local modes of filament deformation, since a previous study [68] in which filaments were represented as connected chains of smaller segments predicted linear dependence of effective viscosity on filament lifetime. While previous studies have emphasized individual roles for cross-link unbinding or filament turnover in stress relaxation [7,13,22], here we have capture their distinct contributions within a single self-consistent modeling framework. Our simulations confirm the theoretical prediction [37,39,69] that spatial heterogeneity of motor activity and asymmetric filament compliance are sufficient to support macroscopic contraction of unconstrained networks. However, under isometric conditions, and without filament turnover, our simulations predict that active stress cannot be sustained. On short timescales, motor forces drive local buildup of extensional stress, but on longer timescales, active local filament rearrangements lead, invariably, to a decay in active stress. These rearrangements can lead to macroscopic network tearing and fragmentation, as previously described [40,47]. However, stress decay can also occur when the distribution of network filaments remains more uniform and the network remains globally connected in the sense that every filament overlaps, and can exchange frictional crosslink forces, with many others. Under these conditions, network rearrangements involve a slower rebalancing of extensile and compressive forces on network filaments. Our results suggest that when filaments can slide relative to one another, the motor forces that produce active stress will drive changes either in connectivity, or in the distributions of forces along individual filaments, or both, that inevitably lead to a decrease in active stress. Thus for contractile networks to maintain isometric tension on long timescales, they must either form stable crosslinks to prevent filament rearrangements, or they must continuously turnover network filaments (or active motors) to renew the local potential for production of active stress. Indeed, our simulations predict that filament turnover is sufficient for maintenance of active stress. As in the passive case, they predict biphasic dependence of steady state stress on filament turnover: For short-lived filaments (τr < τm), steady state stress increases linearly with filament lifetime because filaments have more time to build towards peak extensional stress before turning over. For longer-lived filaments (τr > τm), steady state stress decreases monotonically with filament lifetime because local rearrangements decrease the mean contributions of longer lived filaments. These findings imply that for cortical networks that sustain contractile stress under approximately isometric conditions, tuning filament turnover can control the level of active stress, and there will be an optimal turnover rate that maximizes the stress, all other things equal. This may be important, for example in early development, where contractile forces produced by cortical actomyosin networks maintain, or drive slow changes in, cell shape and tissue geometry [7,70]. For cortical networks that undergo steady state flows driven by spatial gradients of motor activity, our simulations predict that the biphasic dependencies of steady state stress and effective viscosity on filament lifetime define multiple regimes of steady state flow, characterized by different dependencies on filament turnover (and other network parameters). In particular, the approximately linear dependencies of steady state stress and effective viscosity on filament lifetime for short-lived filaments define a fast turnover regime in which steady state flow speeds are buffered against variations in filament lifetime, and are predicted to depend in a simple way on motor activity and crosslink resistance. Measurements of F-actin turnover times in cells that undergo cortical flow [32,71–75] suggests that they may indeed operate in this fast turnover regime. For reference values of model parameters, the steady state strain rates predicted for the high turnover regime (≈ 2x10−4/sec) are approximately ten-fold lower than those measured in polarized C. elegans zygotes (≈ 1 − 2x10−3/sec) [6,11,76]. This is reasonable agreement, given uncertainties about these reference values. For example, our reference value for effective crosslink friction ξ is 10-100-fold higher than friction coefficients measured for single crosslinkers and molecular motors in vitro [77,78]. A ten-fold lower value for this parameter would yield a ten-fold increase in the predicted steady state strain rate (12). Interestingly, recent studies in C. elegans embryos suggests that cortical flow speeds are surprisingly insensitive to depletion of factors (ADF/Cofilin) that govern filament turnover [11], again consistent with our model’s predictions. Stronger tests of our model’s predictions will require more systematic analyses of how flow speeds vary with filament and crosslink densities, motor activities, and filament lifetimes.
summary: In this paper, we develop and analyze a minimal model for a 2D network of cross-linked actin filaments and myosin motors, representing the cortical cytoskeleton of eukaryotic cells. We implement coarse-grained representations of force production by myosin motors and stress dissipation through an effective cross-link friction and filament turnover. We use this model to characterize how the sustained production of active stress, and the steady dissipation of elastic stress, depend individually on motor activity, effective cross-link friction and filament turnover. Then we combine these results to gain insights into how microscopic network parameters control steady state flow produced by asymmetric distributions of motor activity. Our results provide a framework for understanding how local modulation of microscopic interactions within contractile networks control macroscopic quantities like active stress and effective viscosity to control cortical deformation and flow at cellular scales.
keywords: cell motility
  actin filaments
  classical mechanics
  fluid mechanics
  drag
  mechanical stress
  viscosity
  molecular motors
  actin motors
  materials science
  network analysis
  damage mechanics
  materials physics
  motor proteins
  physical chemistry
  chemical properties
  contractile proteins
  computer and information sciences
  proteins
  deformation
  chemistry
  fluid dynamics
  continuum mechanics
  physics
  biochemistry
  cytoskeletal proteins
  cell biology
  myosins
  biology and life sciences
  physical sciences
title: Filament turnover tunes both force generation and dissipation to control long-range flows in a model actomyosin cortex