Quantifying Dissipation in Actomyosin Networks
Carlos Floyd Garegin A. Papoian Christopher Jarzynski






1]Carlos Floyd
1,2,3,†]Garegin A. Papoian
1,2,3,4,*]Christopher Jarzynski
[1]Biophysics Program, University of Maryland, College Park, MD 20742 USA
[2]Department of Chemistry and Biochemistry, University of Maryland, College Park, MD 20742 USA
[3]Institute for Physical Science and Technology, University of Maryland, College Park, MD 20742 USA

[4]Department of Physics, University of Maryland, College Park, MD 20742 USA
[†]email: gpapoian@umd.edu
[*]email: cjarzyns@umd.edu





Quantifying Dissipation in Actomyosin Networks
    [
    December 30, 2023
==============================================

Quantifying entropy production in various active matter phases will open new avenues for probing self-organization principles in these far-from-equilibrium systems.  It has been hypothesized that the dissipation of free energy by active matter systems may be optimized to produce highly dissipative dynamical states, hence, leading to spontaneous emergence of more ordered states. This interesting idea has not been widely tested. In particular, it is not clear whether emergent states of actomyosin networks, which represent a salient example of biological active matter, self-organize following the principle of dissipation optimization. In order to start addressing this question using detailed computational modeling, we rely on the MEDYAN simulation platform, which allows simulating active matter networks from fundamental molecular principles. We have extended the capabilities of MEDYAN to allow quantification of the rates of dissipation resulting from chemical reactions and relaxation of mechanical stresses during simulation trajectories. This is done by computing precise changes in Gibbs free energy accompanying chemical reactions using a novel formula, and through detailed calculations of instantaneous values of the system's mechanical energy. We validate our approach with a mean-field model that estimates the rates of dissipation from filament treadmilling. Applying this methodology to the self-organization of small disordered actomyosin networks, we find that compact and highly cross-linked networks tend to allow more efficient transduction of chemical free energy into mechanical energy. In these simple systems, we do not observe that spontaneous network reorganizations lead to increases in the total dissipation rate as predicted by the dissipation-driven adaptation hypothesis mentioned above. However, whether such a principle operates in more general, more complex cytoskeletal networks remains to be investigated.

 










§ INTRODUCTION



The actin-based cytoskeleton is a dynamic supramolecular structure that, by sustaining and releasing mechanical stress in response to various physiological cues, mediates the exertion of force by cells both on their environments and within their bodies <cit.>.  These cytoskeletal structures are composed of long (on the order of 1 μ m in vivo <cit.>) actin polymers which are interconnected by various cross-linkers, as well as by myosin motor filaments, resulting in a three-dimensional network-like organization referred to as an “actomyosin network" <cit.>.[In our terminology, we will distinguish between “cross-linking proteins," which will include both active (e.g. myosin filaments) and passive (e.g. α-actinin and fascin) proteins that bind to adjacent actin filaments, and “cross-linkers," which refer exclusively to passive cross-linking proteins.] 
Part of the intricacy of actomyosin network dynamics is due to the mechanosensitive kinetic reaction rates controlling cross-linker and myosin filament unbinding as well as myosin filament walking: at high tension, cross-linkers will unbind more quickly (slip-bond) whereas motors will unbind and walk less quickly (catch-bond and stalling). <cit.>.  These reactions control the actomyosin network connectivity, which in turn determines the ability of the network to globally distribute stress <cit.>.  Thus the mechanosensitive feedback introduces nonlinear coupling between the stress sustained by an actomyosin network and the network's ability to reorganize in response to that stress.  In order to be responsive to physiological cues, the dynamics of these systems occur far from thermodynamic equilibrium; the hydrolysis of an out-of-equilibrium concentration of ATP molecules fuels a) the stress-generating activity of the myosin motor filaments, and b) filament treadmilling <cit.>.  Filament treadmilling is a steady-state situation in which the polymerization at the plus end of the filament is compensated by the depolymerization at the minus end, resulting in the filament moving forward without its length changing.  As a result of these these local free energy-consuming processes, actomyosin networks constitute an interesting and biologically important example of soft active matter. Active matter is composed of agents that individually transduce free energy from some external source, in this case the chemical potential energy of many ATP molecules <cit.>.  Dissipation in these systems results when the free energy consumed Δ G is greater than the quantity of work W done by the system on its environment, with the remainder Δ G - W serving to increase the total entropy.         


The viewpoint of actomyosin networks as active matter systems has been fruitfully adopted in recent theoretical and experimental studies, yet a lack of ability to quantify the rates of free energy transduction by these systems has hindered development of some of these lines of study.  The emergence of distinctive dynamical states (for instance pulsing actin waves or vortices) during the self-organization of actomyosin systems has been documented in several in vitro experiments <cit.>.  These emergent patterns depend sensitively on the concentrations of myosin filaments and cross-linkers:  myosin filament concentration controls the level of active stress generation, and cross-linker concentration controls the degree of mechanical coupling of actin filaments, which has been described using the language of percolation theory <cit.>.  While these emergent dynamic patterns have been characterized in detail, a general mechanism explaining why these patterns emerge under given conditions has not yet been proven.  It might be expected, given that these systems operate away from thermodynamic equilibrium, that the quantity of free energy dissipated during a system's evolution is optimized, similar to the principle of minimum entropy production in the near-equilibrium theory of irreversible thermodynamics <cit.>.  However, this minimum entropy production principle breaks down in the far-from-equilibrium, nonlinear-response regime, where many active matter systems including actomyosin networks operate <cit.>.  It has recently been proposed that another optimization principle applies arbitrarily far from equilibrium.  This idea, referred to as dissipation-driven adaptation, suggests that, in general, a coarse-grained trajectory of some nonequilibrium system will be more likely than all alternative trajectories if the amount of free energy absorbed and dissipated along that trajectory is maximal <cit.>.  This organizing principle has been borne out in model systems <cit.>, yet has also been shown to have certain counter examples <cit.>; it remains actively debated.  In the case of actomyosin systems, it has not yet been tested because of the difficulty in measuring dissipated free energy using most experimental approaches.  In this paper we take the first steps toward such a test, by developing a simulation methodology allowing the quantification of dissipated free energy during the self-organization of actomyosin networks.

In addition to being of interest in the field of active matter systems, dissipation in actomyosin networks has also been an important factor in recent experimental studies of cell mechanics.  Rheological properties of actomyosin networks largely determine rheological properties of the whole cell, and it has been discovered that the dissipative component of the cell's viscoelastic response to mechanical oscillations (called the loss modulus) is partly attributable to the ATPase activity of myosin motor filaments in actomyosin cytoskeleton <cit.>.  In the context of cell mechanics, myosin filaments have several roles: they produce mechanical stress by pulling on actin filaments, they dissipate mechanical stress by disassembling actin filaments and higher order stress-sustaining filament structures, and they consume chemical energy through the hydrolysis of ATP molecules.  It has been proposed that a lack of detailed understanding of the effects of these processes, and of dissipation in actomyosin systems more generally, underlies inconsistent, widely variable traction force microscopy measurements of cell migration <cit.>.  Progress along this line is hindered by an absence of methods to study dissipation in actomyosin networks directly and at sufficiently high spatio-temporal resolution.  

To address these needs, we introduce a computational approach to measure dissipation during simulations of actomyosin network self-organization using the simulation platform MEDYAN (Mechanochemical Dynamics of Active Networks) <cit.>.  MEDYAN simulations marry stochastic reaction-diffusion chemistry algorithms with detailed mechanical models of actin filaments, cross-linkers, myosin motor filaments, and other associated proteins, and it also accounts for mechanosensitive reaction rates.  This combination of simulation features makes this software uniquely capable of probing the complexity of actomyosin network dynamics.  For instance, past studies utilizing MEDYAN have investigated the dependence of network collapse on myosin filament and cross-linker concentrations, as well as the origin of local contractility in actomyosin networks <cit.>.  We refer the reader to the paper describing MEDYAN for a detailed discussion of the various aspects of the simulation platform <cit.>, while here we describe an extension of that platform that allows for calculation of the energetics of the chemical and mechanical events occurring during simulation.  We utilize these new capabilities to characterize the dissipation resulting from filament treadmilling, for which we further introduce a mean-field model, as well as from myosin filament walking.  We study both the time-dependence and the distributions of dissipation rates as concentrations of cross-linkers and myosin filaments are varied, observing that transduction of chemical energy to stored mechanical energy is more efficient at denser network organizations.  For these simulations, we first explore systems with “plain" myosin filaments and cross-linkers that are not mechanosensitive, in order to simplify the overall dynamics.  We then introduce their mechanochemical coupling to understand its effect on the observed trends.  We end by discussing how this new methodology can provide a valuable technique to advance the studies of actomyosin networks mentioned above.   




§ METHODS



 §.§ Measuring Dissipation in MEDYAN

We first give a brief overview of the MEDYAN simulation platform.  MEDYAN employs a stochastic chemical evolution algorithm in conjunction with mechanical representations of polymers and cross-linking proteins to simulate the dynamics of networks with active components, including but not limited to actomyosin networks.  The simulation space comprises a grid of reaction-diffusion compartments, inside which chemical species (e.g. unpolymerized subunits or cross-linking proteins) are assumed to be homogeneously distributed without specified locations, and which participate in reactions (e.g. (de)polymerization or (un)binding) according to mass-action kinetics; the species can additionally jump between compartments in diffusion events.  When an unpolymerized subunit polymerizes to or nucleates a filament, it becomes part of the mechanical subsystem, gaining location coordinates in the simulation volume and becoming subject to mechanical potentials depending on its interaction with other mechanical elements.  Through chemical reactions such as myosin filament binding and walking, the mechanical energy of the system changes, and the new net forces are then periodically relaxed in a mechanical equilibration phase, using conjugate gradient minimization.  We fill in salient details of the above overview as they become relevant below.  A user provides input data including system size and simulation length, mechanical parameters (e.g. stretching and bending constants and excluded volume cutoff distances), size of the polymer subunits, energy minimization algorithm parameters, chemical simulation algorithm parameters, choices for the modeling of force-sensitive reaction rates, initial conditions of the filaments (either specified or randomly generated), a list of reacting species and their associated parameters, a list of reactions involving those species and their associated parameters, and a list of desired output information.  The output of a simulation is a set of trajectory files containing information at each time point, which can include positions of the mechanical network elements, tensions on the elements, and copy numbers of the chemical species, among other things.  In Supplementary Information we discuss parameterization of the simulations analyzed in this paper.  MEDYAN is extensible in that it is possible to implement new types of outputs, depending on the experimental needs; in this paper we describe a novel output that reports the changes in the Gibbs free energy of the system.


As a MEDYAN simulation progresses, the Gibbs free energy of the system continually changes due to occurrences of chemical reactions and structural rearrangements of the polymer network.  These processes are driven by an out-of-equilibrium concentration of ATP which fuels filament treadmilling and myosin filament walking.  Dissipation measurement in MEDYAN works by calculating running totals of the chemical and mechanical energy changes. The running totals can then be converted into instantaneous rates by taking the numerical derivative at each time point using the forward difference quotient.  The algorithm for tracking these energy changes is compatible with the following sequence of consecutive procedures that make up one iterative cycle of a MEDYAN simulation <cit.>:

	
	
  * Evolve system with stochastic chemical simulation for time t_min.  
	
  * Calculate the changes in the mechanical energy resulting from the reactions in Step 1.
	
  * Mechanically equilibrate the network based on the new stresses calculated in Step 2.  
	
  * Update the reaction rates of force-sensitive reactions based on the new forces.
      

Dissipation tracking is done by calculating for each of these four steps a change in free energy, and then adding these free energy changes to determine the total change in free energy resulting from each cycle, Δ G_dissipated.  Since, at least in this study, the actomyosin network is not mechanically coupled to any work reservoir external to the simulation volume (i.e. W=0), Δ G_dissipated is indeed dissipated energy <cit.>.  This methodology could be straightforwardly extended to account for work exchanged with an external system in future studies, however.  Step 4 in the MEDYAN simulation cycle does not result in a change in free energy, as explained in Supplementary Information.  The notation chosen for these free energy changes, as well as the direction in which the energy is changing during each procedure, is illustrated in Figure <ref>.  




From this picture, we have the following relations:

    Δ G_dissipated   = Δ G_chem + Δ G_mech
    Δ G_dissipated   = Δ G_chem, dissipated + Δ G_mech, dissipated
    Δ G_chem, dissipated   = Δ G_chem + Δ G_stress
    Δ G_mech, dissipated   = Δ G_mech - Δ G_stress

For the depicted relative position of energy levels, the sign convention is such that all values of Δ G except for Δ G_stress will be negative (indicated by the arrow's direction), since G refers to the free energy of the system, not of its environment, and will therefore be tend to be negative as the system moves down the free energy landscape.  The usual intuition that the total dissipation is positive can be stated 

    Δ G^+_dissipated = - Δ G_dissipated >0 ,

where the superscript “+" indicates the positive change in the total entropy.

Equation <ref> says that, given some change in the system's chemical potential energy, Δ G_chem, resulting from reactions occurring during Step 1, a portion of that energy is used to deform the polymer network (e.g. via myosin filament pulling on actin filaments).  This increases the mechanical energy in the network by an amount Δ G_stress.  Only the portion of Δ G_chem which has not gone into Δ G_stress has been dissipated as heat.  In Step 3 the network is mechanically equilibrated, resulting in relaxation of net forces (though not of all stresses) and updating of the network elements' positions.  We refer to the decrease in mechanical energy resulting from this relaxation as Δ G_mech, dissipated.  

The calculation of Δ G_stress and Δ G_mech, dissipated is based on a set of mechanical potentials describing interactions between elements of the actomyosin network.  Polymers are modeled as a sequence of thin, unbendable, yet extensible cylinders that are joined at their ends by beads whose positions define the polymer's configuration.  The structural resolution of MEDYAN is at the level of the cylinders, which in this study are 27 nm long and have effective diameters of approximately 5 nm, however the diameter is not a parameter of the simulation, being instead effectively determined by the strength of the excluded volume interaction between cylinders.  Cross-linking proteins (e.g. α-actinin and myosin filaments) are modeled as Hookean springs connecting these cylinders by attaching to discrete binding sites.  Included among the mechanical potentials are various modes of filament deformation, excluded volume interactions, and stretching of cross-linkers and myosin filaments.  Mechanical equilibration is accomplished by constrained minimization of the mechanical energy with respect to the positions of the network elements.   A full description of the mechanical potentials and equilibration protocols is given in <cit.>.  Determining Δ G_stress and Δ G_mech, dissipated requires evaluating the instantaneous total mechanical energy of the system at certain points during the iterative simulation cycle and taking the difference of those values.  

The calculation of Δ G_chem and Δ G_chem, dissipated is accomplished by incrementing a running total of the chemical energy G_chem whenever a reaction stochastically occurs during Step 1, and finding the accumulated change at the end of the protocol.  Chemical stochastic simulation in MEDYAN uses a Gillespie-like reaction-diffusion algorithm over a grid of compartments which constitutes the simulation volume.  Diffusing species are assumed to be homogeneous (i.e. obey mass-action kinetics) inside the compartments, and can jump between the compartments leading to concentration gradients at the scale of the compartment length (taken to be roughly the Kuramoto length of diffusing G-actin, following <cit.>).  The evolving polymer network is overlaid on this compartment grid, with each piece of a polymer reacting with diffusing species according to the concentrations in its local compartment.  Again, we refer the reader to <cit.> for a more detailed description of the chemical dynamics.  For the present purpose of measuring dissipation, we introduced into this simulation protocol a precise formula for the change in Gibbs free energy corresponding to the occurrence of various reactions as a function of the instantaneous compartment concentrations.  In the Supplementary Information we establish this formula, while we present its derivation in an accompanying paper <cit.>.  The set of chemical reactions used to describe actomyosin networks in this study is based on a previous model of actin polymerization dynamics that explicitly treats hydrolysis states of the nucleotide bound to each actin subunit <cit.>.  This level of detail allows to quantify the dissipation resulting from ATP hydrolysis during filament treadmilling.  To increase computational efficiency, we neglect the dynamics of nucleotide hydrolysis states of the tips of the filaments.  This has been shown in previous work to be a valid approximation to the full dynamics which includes the states of the filament tips <cit.>.  We refer to the resulting set of reactions describing actin polymerization dynamics as the Constant Tip (CT) model.  However unlike in the original CT model, here we explicitly include G-actin bound to ADP-Pi as a reacting species, for completeness and since the extra computational strain of doing so is small.  The actin subunit species tracked in this model are distinguished by their polymerization state and by the hydrolysis states of the nucleotide to which they are bound.  We notate a species as G or F, to represent globular (un-polymerized) or filamentous (polymerized) actin respectively, superscripted by T, Pi, or D to represent that it is bound to ATP, ADP-Pi, or ADP, respectively; thus for instance filamentous actin bound to ADP-Pi is notated F^Pi. We also include reactions describing cross-linker (un)binding and myosin filament (un)binding and walking.  We exclude filament nucleation, severing, destruction, and annealing reactions, thus the number of filaments is constant throughout the simulation trajectories.  In the Supplementary Information we describe how we compute the change in Gibbs free energy for each reaction in this set, as well as how we parameterize the simulations.

Lastly, we developed a mean-field model to describe just the dissipation resulting from reactions in the CT model, i.e. excluding cross-linkers and myosin filaments.  We describe the model and present its results in Supplementary Information.  We find, among other things, that the steady-state dissipation rate from filament treadmilling counter-intuitively does not depend on the total amount of actin, but only on the number of filaments.   






§ RESULTS





 §.§ Total Dissipation Rates of Disordered Networks Do Not Increase


We studied dissipation rates accompanying the process of myosin-driven network self-organization.  We first excluded in these simulations the force-sensitivity of the reaction rates describing cross-linker unbinding and myosin filament unbinding and walking.  This allowed us to understand a simplified version of the dynamics (we later discuss the the effect of including mechanochemical feedback).  We analyzed the trajectories of the quantities Δ G_dissipated, Δ G_chem, dissipated, Δ G_mech, dissipated, Δ G_chem, and Δ G_mech over a set of simulations with identical initial concentrations but different random filament distributions.  In Figure <ref>, we display these trajectories averaged over 10 runs.  We used a simulation volume of 1 μ m^3, divided into 8 compartments, and initial conditions of equal amounts (10 μ M each) of G^T and G^D actin in a 0.08 μ M pool of seed filaments containing F^T,  as well as 0.2 μ M myosin and 1.0 μ M cross-linkers.  Similar trajectories for other conditions are shown in SI Figures [Due to file size constraints on preprint upload, this supplementary figure has been excluded.  It will be available in the print version coming soon in Interface Focus.].  Points lying outside 3 median absolute deviations (MAD) have been excluded for this visualization <cit.> [For certain statistical aggregations in this paper we use the MAD since it is the estimate of scale most robust against outliers, having a breakdown point of 50%.  However, we use the unscaled version of the MAD and, as a result, do not claim that this statistic is a consistent estimator of the standard deviation <cit.>.  The distributions of most of the quantities of interest are too pathological to allow for a straightforward choice of scaling factor, so we simply us the unscaled version, defined as MAD = median_i  |x_i - median_jx_j |, where x_i are elements in the data set.].  




	



For each tracked quantity, there is an initial transient phase (lasting a few tens of seconds) followed by fluctuations around a roughly steady value; we do not observe in this set of simulations any slow approach to a significantly different total dissipation rate, which might correspond to large reorganizations of the network.  However for high concentrations of cross-linkers (C_CL = 5.0 μ M) and myosin (C_M = 0.4 μ M), we observe that the contribution of Δ G^+_mech, dissipated to Δ G^+_dissipated tends to increase relative to that of Δ G^+_chem, dissipated (SI Figure dns).  It is known that network percolation and collapse occurs only under certain conditions of myosin and cross-linker concentrations <cit.>, and it is under these conditions which result in collapse that we observe an increase in Δ G^+_mech, dissipated relative to Δ G^+_chem, dissipated (SI Figures dns).  This indicates that more mechanical stress is being created by myosin filament walking as the network collapses and becomes more densely cross-linked. The transient phase corresponds to the initial polymerization of the seed filaments followed by the initial mechanical coupling of filaments by cross-linkers and myosin filaments.  Following the transient phase, the networks in this study are generally disordered (Figure <ref>).  The dissipation rate corresponding to the initial polymerization of the seed filaments is much larger than the chemical dissipation resulting from myosin filament activity.  However, following the transient phase, after which the treadmilling dissipation has reached its steady-state (SI Figure <ref>), the contribution from myosin filament activity outweighs the dissipation resulting from filament treadmilling.  Tracking the instantaneous rate of change in Δ G^+_chem resulting from each reaction separately, we observe myosin filament walking to contribute the majority to Δ G^+_chem after the transient phase, however the amount depends on C_CL and C_M.  The integrated contributions of each reaction to Δ G^+_chem for different conditions are shown in Figure <ref> and SI Figures dns.  Interestingly, diffusion contributes an appreciable fraction to Δ G^+_chem; plus and minus ends of actin filaments tend to localize together (Figure <ref>) and through treadmilling deplete the concentrations of G^T and G^D in certain reaction compartments relative to others, leading to significant diffusion gradients on the scale of the compartment length.  This suggests that the establishment of concentration gradients is an important driving force in actomyosin self-organization, as has been noted in other work <cit.>.  In the initial transient phase, the mechanical energy changes appreciably as the filaments grow and are initially coupled to each other.  Following this, however, the rate of Δ G^+_mech is on average near zero.  This indicates that, despite the process of mechanical stress generation through myosin filament activity and treadmilling, the resulting stress is dissipated through fast relaxation such that, on a slower timescale, the mechanical energy of the system does not change in a significant, persistent way.     



	





	


 
Aggregating each of the 2000 s of the 10 trajectories into a collective data set for each condition of C_M and C_CL, we next analyzed the distribution of the instantaneous rates of Δ G^+_dissipated, Δ G^+_mech, dissipated, Δ G^+_chem, dissipated, Δ G^+_stress, Δ G^+_mech, and Δ G^+_chem.  In Figure <ref> we plot histograms of these 6 quantities for the conditions C_CL = 1.0 μ M, C_CL = 0.2 μ M.  In SI Figures dns we show the same plots for other conditions.   



	



Each distribution contained heavy tails which, for the purpose of fitting, we suppressed by excluding data lying outside 5 MAD's of the median.  No distributions were sufficient to cleanly fit the full set of data, so we focus here on the center of the distribution and include a qualitative discussion of the heavy tails below. With the exception of Δ G^+_mech which was fit with a generalized normal distribution, each distribution exhibited significant skew and could be fit reasonably well with a log-normal or a gamma distribution.  At higher concentrations of myosin and cross-linkers, the histograms were less cleanly fit by any standard distributions (SI Figure dns).


Log-normal distributions are fairly ubiquitous across different fields and systems <cit.>.  For instance, it has been shown that the distributions of concentrations of species in a chemical reaction network are log-normal  <cit.>.  Gamma distributions are similarly common and often difficult to discriminate from log-normal distributions <cit.>.  A speculative explanation for the gamma distribution of Δ G_stress is as follows: Δ G_stress can be viewed as resulting from a number of myosin filament steps that, according to the central limit theorem, is approximately normally distributed given a sufficiently long time between simulation snapshots t_snap (these stepping events are not truly i.i.d., but to a first approximation we may assume they are).  The main effect of each of these steps is to increase the harmonic stretching potential on the myosin filaments as well as on the actin filament cylinders by a roughly fixed stepping distance.  Thus the increase in mechanical energy, Δ G_stress, is approximately a quadratic function of the normally distributed number of myosin filament steps per t_snap.  As shown in SI Figure dns, the resulting distribution of Δ G_stress is well-fit by a gamma distribution and bears qualitative similarity to the histogram of Δ G_stress in Figure <ref>.  In Figure <ref>, Δ G_stress and Δ G^+_mech, dissipated have similar distributions, indicating that, for each MEDYAN cycle, almost all the stress accumulated following Step 1 is then immediately relaxed.  The remainder goes into Δ G^+_mech, whose distribution is centered on zero with little skew.  Furthermore,  the distribution of Δ G^+_mech, dissipated has particularly heavy tails, indicating infrequent yet large relaxation events.  This tendency has some precedent in avalanche-prone systems, whose hallmarks are self-organized criticality, intermittency, and scale-invariance in their distribution of avalanche event sizes <cit.>.  While we do not observe true scale-invariance in this set of simulations (i.e. power-laws cannot fit these distributions cleanly), it will be interesting to continue to explore actomyosin networks dynamics in the framework of self-organized criticality <cit.>.






 §.§ More Compact Networks are More Efficient

We simultaneously varied concentrations of myosin filaments C_M, and cross-linkers, C_CL, which is known to produce a range of network architectures <cit.>.  Using cross-linker concentrations of 0.1, 1.0, and 5.0 μ M, and myosin concentrations of 0.1, 0.2, and 0.4 μ M, we studied the effects on Δ G^+_dissipated, Δ G^+_chem, dissipated, and Δ G^+_mech, dissipated.  The concentrations of actin subunits and filaments are the same as described above.  For each of the 9 conditions, we ran 10 simulations of 2000 s.  In Figure <ref>, we display the median values of Δ G^+ _chem, dissipated and Δ G^+ _mech, dissipated along each trajectory and over each repeated trajectory for these conditions.  The median is used here because of its insensitivity to outliers such as are found in the heavy-tailed distributions of these quantities <cit.>.  
  


	



We find that as C_M is increased, the median rates of Δ G^+_dissipated, Δ G^+_chem, dissipated, and Δ G^+_mech, dissipated all tend to increase.  Further, the value of of Δ G^+_mech, dissipated relative to Δ G^+_chem, dissipated increases.  As C_CL is increased with C_M fixed, then Δ G^+_chem, dissipated tends to decrease.  Increased concentrations of myosin filaments obviously have a strong effect on the dissipation simply because they are the chief active agents in the system once treadmilling has reached its steady-state.  We can define a measure of efficiency for the present purpose as:


    η = Δ G_stress/Δ G^+_chem = Δ G^+_mech, dissipated - Δ G^+_mech/Δ G^+_dissipated - Δ G^+_mech≈Δ G^+_mech, dissipated/Δ G^+_dissipated

where the approximation follows since, as illustrated in Figure <ref>, the average of Δ G^+_mech is close to zero.  Thus it is evident that, perhaps surprisingly, as C_M increases, η increases: the more motors are in our system, the more efficiently can chemical energy be converted into mechanical stresses.  Further, as C_CL increases, η tends to increase because Δ G_chem, dissipated decreases relative to Δ G_mech, dissipated.  At higher levels of cross-linking, which introduce mechanical constraints on the filaments, the walking of myosin motor filaments will produce more stress compared to at low levels of cross-linking, when filament sliding can result from myosin filament walking, producing less stress.  

We repeated the above experiments with the inclusion of mechanochemical feedback on the reaction rates controlling cross-linker unbinding and myosin filament unbinding and walking  (SI Figures dns).  We observed, somewhat surprisingly, no qualitative differences compared to the results described for the case of no feedback, however there were quantitative differences in the stress production and radius of gyration for certain conditions.  These quantitative differences result from the fact that, for the concentrations C_M and C_CL used above, we did not observe significant collapse of the actomyosin network when mechanochemical feedback was included due to the stalling and catching of the myosin filaments.  As a result the networks were less densely cross-linked, and the efficiency was lower.  The total dissipation rates are largely unaffected by the inclusion of mechanochemical feedback; it was primarily the degree to which chemical energy had been converted to stress that was different for certain conditions.  


It is worthwhile to mention how these results compare with in vitro studies of dissipation in actomyosin systems.  While the computational approach described in this paper allows uniquely highly-resolved and direct measurement of free energy changes, other experimental methods have produced qualitatively similar results to those obtained here.  Rheological experiments have determined that a single cell's response to compression is similar in nature to a muscle's response to increasing load, suggesting that the actomyosin network underlies the cell's mechanical responsiveness <cit.>.  Further, this responsiveness is modulated by blebbistatin, a myosin ATPase inhibitor, highlighting myosin's role in negotiating how the network rearranges in response to the sustained stress.  Using the metric of mechanical dissipation to measure the degree of structural rearrangements and release of stress, we confirm that these processes indeed sensitively depend on myosin activity, which we control here through its concentration.  Additional rheological studies have probed the mechanical dissipation of actin cortices more directly, using the loss modulus as a readout.  These have also indicated that inhibiting myosin activity reduces the mechanical dissipation of the system, causing it to be behave more elastically <cit.>.  Lastly, we mention a recent study that quantified mechanical dissipation of actin filaments using a novel experimental method <cit.>.  By measuring the flow through a low-dimensional phase space defined by the amplitudes of the filament bending modes <cit.>, they determine the entropy production of fluctuating actin filaments in different phases of contractility.  They relate the entropy production to the degree of transverse bending of the filaments, as opposed to sarcomeric filament sliding, caused by myosin filament walking.  Similarly, we here relate the hindrance of filament sliding due to cross-linker density to increased mechanical dissipation rates.  Quantitative comparisons of these two approaches to mechanical dissipation measurement will be an interesting future direction.  We note finally that a unique capability of quantifying dissipation using MEDYAN is the ability to simultaneously measure the energetics of chemical reactions in addition to the changes in the mechanical energy, which is not currently available using in vitro methods.  










§ DISCUSSION


We have introduced a methodology for tracking the energetics of  chemical and mechanical events during a MEDYAN simulation, allowing us to probe the properties of actomyosin networks as dissipative active matter systems.  The distinction between dissipation's mechanical and chemical origins is natural in the context of MEDYAN's iterative simulation procedure which carries out chemical stochastic simulation, mechanical deformation, and mechanical relaxation at separate times.  As explained in <cit.>, this procedure exploits a separation of timescales between the characteristic mechanical relaxation times of actomyosin networks <cit.> compared to typical waiting time between reactions that introduce mechanical stresses, such as myosin filament walking <cit.> or filament growth <cit.>.  Ultimately the source of all dissipation is the chemical potential of ATP molecules driving treadmilling and myosin filament activity.  This is reflected by the near equality of Δ G^+_chem and Δ G^+_dissipated in Figure <ref> and SI Figures dns.  On a fast timescale (that of the characteristic mechanical relaxation time), however, the free energy of chemical reactions cause small force deformations of the network which are then quickly and almost fully relaxed.  Thus, for a myosin filament stepping event, only a portion of the chemical free energy Δ G^+_chem, dissipated is immediately dissipated as heat, with the rest going into temporarily increased mechanical energy of the actomyosin network, Δ G_stress.  The fast relaxation of this new mechanical stress constitutes what we refer to as mechanical dissipation, Δ G^+_mech, dissipated, and the small residual stress after this relaxation has balanced all net forces acting on the system results in a change of the mechanical energy of the system on a slow timescale, Δ G^+_mech.
 

One interpretive framework which is useful to understand the flow of free energy in actomyosin networks is illustrated in Figure <ref>.  We can think of the different forms of free energy storage, including chemical potentials, mechanical stress, concentration gradients (which could also be considered as arising from chemical potential differences across compartments), and dissipated energy, as nodes on a directed graph, where edges represent transduction of energy from one form of storage to another.  The weights of these edges represent the amount of free energy flowing through them.  In this picture, we can describe the process of network percolation, which occurs at increasing concentrations of cross-linkers and myosin filaments, as widening the edge flowing from chemical potential into mechanical stress, while thinning the edge from chemical potential directly to dissipation.  The edge weights corresponding to the establishment of concentration gradients and the resulting diffusive dissipation will not be affected dramatically by the onset of percolation except to the extent that percolated networks might lead to the formation of more bundles, and therefore more significant concentration gradients.  At steady state, we have a stationary current on the graph, fueled by the chemical potential of the assumed limitless supply of ATP.  Of course, at this level of description, we have coarse-grained away the details of the specific chemical reaction networks and mechanical potentials which constitute the system, however by doing so we gain a clearer understanding of how percolation of the network alters the flows of free energy.   



	




The new capabilities of MEDYAN allowing detailed energetics computations should provide a way to address outstanding issues in the the fields of active matter systems and of cell mechanics.  In this study we observe dissipation rates to stay at a low fluctuating steady-state following a high initial transient phase, apparently in contrast with the predictions of the dissipation-driven adaptation hypothesis.  However we also do not observe slow reorganization of the actomyosin network into different higher order structures such as bundles (for the concentrations of components used here we observe disordered networks only), and as a result we can not rule out that such reorganizations correspond to marked changes in the dissipation rates.  A dedicated test of the hypothesis of dissipative-driven adaptation should be straightforward with this methodology.  For instance, a simulated gliding assay, which has been shown in vitro to lead to diverse dynamical patterns <cit.>, may indicate that these emergent patterns correspond to an optimization of the chemical dissipation from myosin filament walking, as argued in <cit.>.  In the context of cell migration, studies investigating the mechanically dissipative activity of myosin filaments are also feasible if one incorporates in simulation an external substrate against which the actomyosin network pulls.  In this type of study, it should be straightforward to determine how the amount of stress which is sustained against the substrate is altered by myosin activity given the energetics calculations in the methodology described in this paper.  It is also feasible to do simulated measurements of the dynamic shear modulus by compressing the actomyosin network at different frequencies.  It should then be possible to directly observe the degree to which the dissipation of elastically stored energy is attributable to myosin walking.  In fact, this last research question has already been investigated to some extent using an alternate computational model to the one presented here <cit.>.  We hope some of these potential future experiments will shed light on outstanding questions in the studies of actomyosin networks.  





§ SUPPLEMENTARY INFORMATION FOR “QUANTIFYING DISSIPATION IN ACTOMYOSIN NETWORKS" 





	

 §.§ Δ G of Chemical Reactions

	
	
	For a reaction of the general form 
	
    ν_1 X_1 + ν_2 X_2 + …⇌υ_1 Y_1 + υ_2 Y_2 + …

	we define X_i as the reactant species, Y_i as the product species, and ν_i and υ_j as their stoichiometric coefficients.  The “stoichiometric difference" is defined as
	
    σ = ∑_j ∈ Pυ_j - ∑_i ∈ Rν_i
 
	where P is the set of products and R is the set of reactants.
	We further define the conversion factor between the copy number of species i, N_i, and its concentration C_i,
	
    Θ = N_Av V

	where N_Av is Avogadro's number, V is the volume of the compartment where the reaction occurs, and we have N_i = C_i Θ.  Next we define the following quantity, reminiscent of the reaction quotient:
	
    Q = ∏_i ∈ R(N_i - ν_i)^N_i-ν_i/N_i^N_i∏_j ∈ P(N_j + υ_j)^N_j+υ_j/N_j^N_j

	where R is the set of reactant species and P is the set of product species.  Lastly, defining Δ G^0 as the standard state change in Gibbs free energy, we arrive at the following expression for the change in Gibbs free energy as a function of the instantaneous vector of species copy number N:
	
    Δ G(N) = Δ G^0 - σ k_B T  logΘ - σ k_B T + k_B T logQ.

	Although this equation might appear unusual, upon making some approximations leveraging the relative size of the stoichiometric coefficients and the species copy numbers it can be shown to reduce to the familiar textbook expression
	
    ΔG(C) = Δ G^0 + k_B T logQ

	where 
	
    Q = ∏_i ∈ R C_i^-ν_i∏_j ∈ P C_j^υ_j

	is the reaction quotient of the species concentrations C_i, and C is the vector of these concentrations.  In the context of a compartment-based reaction diffusion scheme, the species copy numbers used in Equation <ref> are specific to the compartment in which the reaction was drawn.  We prefer to use Equation <ref> in simulation because small deviations from this nearly exact expression can actually lead to significant systematic bias of the change in Gibbs free energy resulting from certain reactions.  This is especially true for reactions that are very frequent, such as the diffusion reaction between adjacent compartments.   Δ G for diffusion reactions, in which a molecule jumps from a compartment where its concentration is N_i,A to a compartment where its concentration is N_i,B, is calculated using an expression similar to Equation <ref>: 
	
    Δ G = k_B T log(N_i,A-1)^(N_i,A-1)/N_i,A^N_i,A(N_i,B+1)^(N_i,B+1)/N_i,B^N_i,B.

	Similarly this expression reduces to the familiar formula 
	
    Δ G = k_B T logN_i,B/N_i,A

	upon leveraging the sizes of N_i,A and N_i,B relative to 1.  The derivation of the above results can be found in an accompanying paper <cit.>.
	
	For some of the reactions in the CT model, Equation <ref> can be directly applied and straightforwardly cast into the form of Equation <ref> to give more exact results.  However other reactions including (de)polymerization, myosin filament walking, myosin filament (un)binding, and cross-linker (un)binding, require some additional treatment. 
	
	For (de)polymerization reactions, care should be taken in defining the concentrations of the reactants and products since those molecules include heteropolymers; we would like to avoid requiring that a specific sequence of distinct subunits constitutes a unique chemical species.  This is possible to do because the polymerization of actin subunits is independent of the chemical identity of the subunit at the tips of the polymer, and is therefore independent of the chemical identity of any of the polymerized subunits <cit.>.  Then, polymerization potentially only depends on the polymer length, such that a general reaction for reversible polymerization can be written as 
	
    (m)_n + m ⟷ (m)_n+1

	where (m)_n is a polymer with degree of polymerization n and m is a subunit.  For this reaction Equation <ref> reads
	
    Δ G = Δ G^0 + kT logC_(m)_n+1/C_(m)_nC_m.
  
	Following <cit.>, we make the simplifying assumption that a polymer's chemical reactivity does not depend on its degree of polymerization, provided that the polymer is sufficiently long and cooperative effects do not apply.  With this, Equation <ref> reduces to 
	
    Δ G = Δ G^0 + kT log1/C_m.
  
	This assumption agrees more or less with intuition: the polymerizing subunit does not really “see" the degree of polymerization of the polymer, which is reflected by constant rates of polymerization for polymers of varying lengths. 
	
	Cross-linking proteins are incorporated in MEDYAN as diffusing species that bind and mechanosensitively unbind to pairs of actin filaments, mechanically coupling them.  We treat the change in free energy of the (un)binding reactions analogously to the (de)polymerization reactions.  Ignoring the chemical identities of the heteropolymer subunits, a general cross-linking protein (un)binding reaction can be written as 
	
    (m)_i + (m)_j + L ⟷ (m)_i L (m)_j

	where L is a cross-linking protein and (m)_i is a polymer with degree of polymerization i.  According to Equation <ref>, the change in Gibbs free energy upon this reaction occurring to the right is
	
    Δ G = Δ G^0 + k T logC_(m)_iL(m)_j/C_(m)_iC_(m)_jC_L.

	Similarly to our assumption that the chemical reactivity of polymers is independent of degree of polymerization for sufficiently long polymers, we assume here that the binding affinity of cross-linking proteins is independent of the degree of polymerization and of the number of cross-linking proteins already bound to the pair of filaments, allowing us to simplify Equation <ref> to 
	
    Δ G = Δ G^0 + k T log1/C_L.

	
	The kinetics of myosin filament walking in MEDYAN are based on the Parallel Cluster Model (PCM) <cit.>.  Taking the rate constants of individual myosin head (un)binding reactions as inputs, and accounting for the statistical distribution of the number of bound heads as well as dependence on the force exerted on the filament, the results of the PCM allow us to write kinetic parameters describing the entire myosin filament, including filament (un)binding and walking rates.  In the MEDYAN implementation, each step of a myosin filament can represent several steps of the constituent myosin heads, where each head step represents the completion of a single myosin head cross-bridge cycle <cit.>.  The head steps have a fixed length, d_step, set by experimental measurements of the myosin isoform of interest.  The length of the filament steps, d_total, is determined from the following MEDYAN parameters: the equilibrium length L_cyl of the cylinders that comprise the coarse-grained representation of actin filaments, and the number of binding sites per cylinder N_bs, giving d_total = L_cyl/N_bs.  The binding sites represent the discrete locations on the cylinders which can be occupied by cross-linkers or myosin filaments.  To account for the discrepancy between d_step and d_total, we multiply the filament walking rate by the ratio
	
    s = d_step/d_total.

	When a filament step occurs in MEDYAN, it thus represents the completion of s^-1 myosin head cross-bridge cycles, each of which has the effect of converting one solvated ATP molecule into solvated Pi and ADP.  The Gibbs free energy the filament walking reaction is then
	
    Δ G = s^-1( Δ G^0 + k_B T logC_ADPC_Pi/C_ATP).

	where Δ G^0 refers to the standard change in Gibbs free energy for hydrolysis of ATP.  This expression for Δ G could be further multiplied by a parameter ζ representing the coupling of the ATP hydrolysis cycle to the forward step of the myosin head; here we follow the assumption of tight coupling, i.e. ζ = 1  <cit.>.  
	
	The mechanochemical updating of force-sensitive reaction rates k_± alters the equilibrium constant via K_eq = k_-/k_+, and therefore the change in free energy Δ G corresponding to that reaction via Δ G^0 = k_B T log K_eq.  For a reaction whose rates have been updated due to some applied force F, it can be shown that the new value of Δ G^0 is approximately given by
	
    Δ G^0(F) = Δ G^0|_F=0 + g(F)

	where G^0|_F=0 is the original (i.e. zero-force) value of Δ G^0 and g(F) is equal to the increase mechanical energy due to the applied force <cit.>.  In this modeling, the extra energy g(F) is counted as a part of Δ G_mech, dissipated, not Δ G_chem, dissipated.  So when, for example, a cross-linking protein unbinds under tension F, the zero-force value Δ G^0|_F=0 is used when computing the change in free energy for that reaction, and when mechanical equilibration next occurs, the released stretching energy g(F) is included in the calculation of Δ G_mech, dissipated. 
	
	

 §.§ Parameterization

	Parameterization of the CT model for the purpose of tracking free energy changes during simulation trajectories consists of choosing values of the rate constants (kinetic parameters) and of Δ G^0 (thermodynamic parameters) for all reactions in the model.  Wherever possible, values from the literature are used.  Experimental measurements have determined rate constants for every reaction, however for some reactions the value of Δ G^0 hasn't been reliably measured, to the best of our knowledge.  Below, we describe a technique to solve for these unknown values.
	
	For reversible reactions, where the forward and reverse rate constants k_+ and k_- are known, such as (de)polymerization and (un)binding of cross-linkers, Δ G^0 can be found from 
	
    Δ G^0 = k_B T log K_eq = k_B T logk_-/k_+

	Literature values for the equilibrium constants or of Δ G^0 are used for irreversible reactions, for which k_- is often too small to determine from direct measurement.  Irreversible reactions in this system include all reactions except for (de)polymerization reactions.  For reactions for which literature values of Δ G^0 are unavailable, it is possible to solve for Δ G^0 values based on a self-consistency condition <cit.>: the sum of the Δ G^0 values around a closed loop of reactions in which the number of molecules has not experienced a net change must be zero since the free energy is a state function (equivalently, by Equation <ref>, the product of equilibrium constants around any such loop must be equal to one).  Writing several such closed loops of reactions leads to a system of equations that can be solved for the unknown variables.  Not all possible loops result in independent equations, but we were able to determine the values of two unknown parameters using the loops illustrated in SI Figure <ref>.  
	
	
	
	The loops in SI Figure <ref> imply the following independent system of equations:
	
    K^T K_hyd K_rel K_nex = K^Pi K_ATP
    
    		K^T K_hyd K_phos K_nex = K^D K_ATP

	where K_hyd represents the hydrolysis of ATP by F^T, K_phos represents the release of Pi by F^Pi, K_nex represents nucleotide exchange converting G^D to G^T, K_ATP represents the hydrolysis of ATP in solution producing ADP and Pi, and K_rel represents the release of phosphate by G^Pi.  
	
	Values from the literature <cit.> can be used to determine 6 of the 8 variables in Equations <ref> and <ref>, which thus represent two equations in two unknowns: K_rel and K_nex.  The resulting parameters are listed in Table <ref>.  
	
	Note that it is possible to draw loops such as those in SI Figure <ref> that would imply that the equilibrium constants for polymerization and depolymerization of, for example G^T, should be the same at the plus and minus ends of the filaments.  This condition is not borne out by the experimental values of these equilibrium constants, and this discrepancy is a recognized outstanding problem <cit.>.  Here, we use the literature values for these equilibrium constants and employ the reaction loop method only to determine the parameters K_rel and K_nex.
	
	For reversible binding of myosin filaments, results from the PCM are used to describe binding and unbinding rates, and therefore K_eq  <cit.>.  The filament binding rate is given as 
	
    k_fil, bind = N_t k_head, bind,

	and the unbinding rate is 
	
    k_fil,unbind = N_t k_head, bind[ ( k_head, bind+k_head, unbind/k_head, unbind)^N_t -1 ]^-1,

	where k_head, bind and k_head, unbind describe the binding kinetics of a single myosin head, and N_t is the number of heads in the filament <cit.>.  The resulting expression for Δ G^0 is
	
    Δ G^0 = k T k_fil, unbind/k_fil, bind = - kT log[ ( k_head, bind+k_head, unbind/k_head, unbind)^N_t -1 ].

	
	
	
	We assume that we have chemostatted concentrations of ATP, ADP, and Pi, which we account for implicitly via the effect of these concentrations on the kinetic and thermodynamic parameters of certain reactions.  Thus these species are not explicitly tracked.  The concentrations of these species affect the change in Gibbs free energy associated with the following reactions:
	 
		
  * myosin filament walking
		
  * nucleotide exchange (G^D → G^T)
		
  * phosphate release by F and G-actin (F^Pi→ F^D, G^Pi→ G^D)
	
	The effect of the concentrations of ATP, ADP, and Pi for the these reactions is to simply change the reaction quotient Q which changes Δ G via
	ΔG = Δ G^0 + k_B T logQ.  For instance the nucleotide exchange reaction can be written explicitly as 
	
    ATP + G^D → ADP + G^T
 
	and the change in Gibbs free energy is 
	
    Δ G = Δ G^0 + k_B T log( C_ADP C_G^T/C_ATPC_G^D).

	To treat the concentrations of ATP and ADP implicitly, we rewrite the reaction as 
	
    G^D → G^T
 
	for which the change in Gibbs free energy is 
	
    Δ G = ΔG^0'  + k_B T log(  C_G^T/C_G^D),

	where 
	
    Δ G^0' = Δ G^0 +   k_B T log(  C_ADP/C_ATP).

	A similar approach is taken for the other reactions mentioned above.
	
	
	The concentration of just ATP (since neither ADP or Pi appear implicitly as reactants in any of the reactions of the CT model) affects the kinetics of the following reactions:
	 
		
  * myosin filament walking
		
  * nucleotide exchange
	
	To understand the effect of C_ATP on the myosin filament walking rate, we employ the five state cross-bridge model of a single myosin head described in <cit.>.  In that model unbinding of a head from the filament substrate occurs via two pathways: a slip path, with rate k_35, and a catch path, with an effective rate k_345.  The catch path is a two-step reaction: the release of ADP with rate k_34, followed by the unbinding of the filament head and binding of ATP with rate k_45 = k_T C_ATP.  The effective rate constant of the approximate one-step representation of this reaction is 
	
    k_345 = k_34 k_T C_ATP/k_34 + k_T C_ATP.

	Unless C_ATP is very low, this reaction rate is limited by k_34.  Because the head can unbind by the catch or slip pathway, the rate for unbinding is 
	
    k_head, unbind = k_35 + k_34 k_T C_ATP/k_34 + k_T C_ATP

	The slip path should also be dependent on ATP concentration, since in state 5 the head is ATP-bound, however we follow the authors of <cit.> in neglecting this dependence since the slip path only becomes active under large load.
	
	The rate of the nucleotide exchange reaction also depends on C_ATP and also occurs in two steps.  The full reaction, with ATP and ADP explicitly included, is 
	
    G^D + ATP → G^* + ATP + ADP → G^T + ADP

	where G^* represents actin with no bound nucleotide.  Following <cit.>, this reaction is approximated as a one-step irreversible reaction with ATP and ADP included explicitly, Equation <ref>. We can write the rate k_nex of this approximate reaction as:
	
    k_nex = k_D→*k_*→ TC_ATP/k_D→*+k_*→ TC_ATP≈ k_D→*

	where k_D→* is the dissociation rate of ADP, k_*→ T is the second order rate constant of ATP association, and the approximation holds except at low concentrations of ATP when ADP dissociation is no longer the rate-limiting step.  The values of k_*→ T and k_D→* are presented and discussed in <cit.>.
	
	
	
	
	
	^a - Values of Δ G^0 determined via Equation <ref>.
	
	^b - Δ G^0 determined from K_eq given in <cit.>.
	
	^c - Δ G^0 determined using constraints as described above.
	
	^d - Parameters describing the myosin filament obtained via Equations <ref>, <ref>, <ref>. 
	
	^e - Depends implicitly on C_ATP, C_ADP, and C_Pi; given value applies to standard state.  
	
	^f - Calculated in simulation using results of PCM, see <cit.>.
	
	
	In all the studies in this paper, implicit nucleotide concentrations are taken to be C_ATP = 8 m M, C_ADP = 7 μ M, and C_Pi = 1 m M, corresponding roughly to the amounts found in human muscle after exercise <cit.>.  Other parameters of the system, including mechanical constants, diffusion rates, screening lengths, and boundary cutoffs have been set to the same values given and discussed in <cit.>.  Cylinder equilibrium lengths L_cyl are chosen as 27 nm with 4 binding sites per cylinder for myosin filaments and 1 binding site per cylinder for cross-linkers, giving approximately physiological values for stepping distances of myosin motor filaments and spacing along actin filaments of bound α-actinin.  We note that the form of the mechanochemical models has been changed from those used in <cit.>; the modeling used here is current as of MEDYAN v3.2, and we refer readers to the documentation at http://www.medyan.org/ for further details. 
	
	
	
	
	
	
	

 §.§ Mean-Field Model of Treadmilling Dissipation

	
	To validate the methods for quantifying dissipation using MEDYAN against a simpler representation of actin filament treadmilling, we developed a mean-field description of the dissipation resulting from chemical reactions in the CT model.  Mean-field models of the trajectory of the vector of concentrations of species, C(t), have been formulated previously as an 11-dimensional system of ordinary differential equations (ODEs) (<cit.>), and in the CT model as a 5-dimensional system of ODEs  (<cit.>).  These models describe the polymerization of a concentration N_fil of actin filaments in a pool of actin subunits of total concentration M.  The subunit species tracked by these models are distinguished by their polymerization state and by the hydrolysis states of the nucleotide to which they are bound. The meaning of “mean-field" in this context is the assumption that reacting species are well-mixed over the entire system volume, therefore obeying mass-action kinetics and deterministic dynamics which can be represented by ODEs.  For an instantaneous value of C, we define a function D_Λ(C) representing the instantaneous rate of dissipation due to a set of reactions Λ.  Thus a solution of a mean-field C(t) = N(t)/Θ, allows us to construct the trajectory of the dissipation rate, D_Λ(C(t)).  This function cannot capture the dissipation due to the activity of myosin filaments or cross-linkers, or from relaxation of mechanical stress, because these aspects are not included in the these mean-field models describing filament treadmilling.  The benefit of such a mean-field model is that one can perform systematic variation of parameters with limited computational demands, and we use it here to study the effect of the parameters N_fil and M on the dissipation due to filament treadmilling.
	
	The function D_Λ(C(t)) can be written as a sum over the reactions λ∈Λ of the instantaneous rate of change of the solution's Gibbs free energy due to that reaction:
	
	
    D_Λ(C(t)) = ∑_λ∈ΛΔ G_λ(C(t)) r_λ(C(t)).

	The expressions for Δ G_λ(C(t)) for different reactions are described above.  The instantaneous rate of reaction λ, r_λ(C(t)), is written as usual for mass-action kinetics as
	
    r_λ(C(t)) =Θ k_λ∏_i ∈ R C_i^ν_i

	where R is the set of reactant species for reaction λ, and where we have included the conversion factor Θ to convert r_λ(C(t)) to units of s^-1.
	
	To facilitate the study of dissipation due to filament treadmilling, we define a set of reactions Λ in the CT model which constitutes the dominant cycle a subunit undergoes in the treadmilling process: a) polymerization of G^T to the plus end, b) hydrolysis of ATP by F^T, c) release of Pi by F^Pi, d) depolymerization of G^D from the minus end, and e) nucleotide exchange converting G^D to G^T.  We refer to this set of 5 reactions as the main treadmilling pathway (MTP).  Alternative sequences of reactions whose net effect is similarly the conversion of one molecule of ATP to ADP and Pi are considered as of secondary importance and not included in this analysis.
	
	
	This mean-field description of entropy production can be compared to the description of entropy production rates for chemically reactive systems that emerged from the Brussels school of thermodynamics.  <cit.>.  The results of that school include the minimum entropy-production principle applicable in the linear regime <cit.>, and the general evolution criterion applicable even in the the non-linear regime <cit.>.  Their formalism typically considers the total entropy production rate as a volume integral over the local entropy production rate density, which is itself written as a sum over fluxes multiplied by their corresponding thermodynamics forces defined at each point of the system.  This sum is decomposed into terms representing diffusion and terms representing chemical reactions, and the terms representing the chemical reactions are written such that the fluxes reflect the net reaction rate at that point in the system.  In contrast, our mean-field description neglects concentration gradients and resulting diffusion fluxes, as we assume a homogeneous distribution of the chemical species.  Equation <ref> then represents only the chemical contribution to the entropy production, as a sum over the rates and affinities of the reactions in the system, implicitly integrating over the homogeneous system volume.  We treat the forward and reverse direction for some chemical reaction as separate terms in Equation <ref>, so these rates cannot be considered fluxes which would include the reverse rate as well.  This allows for more general sets of reactions which might include effectively irreversible processes for which the reverse rate is negligible.  However the set of reactions Λ could be chosen to include a reverse reaction for each forward reaction, with the result that these pairs of terms represent fluxes along the reversible reaction pathway.  The parsimonious reaction set MTP is chosen not to fully describe the rate of entropy production in the system, but to allow easy visualization of the main contributions to the entropy production.  We do not pursue the connection of our treatment to the formalism of the Brussels school further here, but we lastly note that our results are compatible with their minimum entropy-production principle, as shown in SI Figure <ref>.  
	
	
	We first verified that the mean-field model of MTP dissipation agreed with results from MEDYAN simulations, to illustrate consistency between these approaches.  In SI Figure <ref> we display the close match between the trajectory of MTP dissipation over a 2000 s run from these two approaches. Note that, to allow direct comparison, only the changes of Gibbs free energy resulting from reactions in the MTP set are visualized for both approaches here, i.e. the contribution from diffusion and other non-MTP reactions in the MEDYAN simulation are not visualized.  We also turned off in MEDYAN the force-sensitive decrease in polymerization rate when the filament tips push against the simulation hard-wall boundaries, since this effect is not captured in the mean-field modeling <cit.>.  We used a simulation volume of 1 μ m^3 and initial conditions of equal amounts (10 μ M each) of G^T and G^D actin in a 0.08 μ M pool of seed filaments containing F^T.  The dissipation rate decreases nearly monotonically, attaining a minimal steady-state value after tens of seconds.  In SI Figure <ref> we also display the individual contributions to the sum in Equation <ref>.  Initial polymerization of G^T to the plus end constitutes the majority of the initial dissipation.  As this polymerization process slows after about 1 second, the hydrolysis of ATP by the now relatively abundant F^T becomes the dominant contribution.  As hydrolysis then slows after about 10 seconds, the total dissipation rate reaches a steady-state value of roughly 80 k_B T/s.  In SI Figure <ref> we also plot the mean-field prediction of the trajectory of the reacting species' concentrations. 
	
	

		
	

	
	We next simultaneously varied the total concentration of actin M and the concentration of actin filaments N_fil, and determined the total dissipation integrated along each trajectory as well as the steady-state dissipation rate.  As M was varied, we held the initial concentration of each species proportionally the same: 49 % G^T, 49 % G^D, and 2 % F^T.   As shown in SI Figure <ref>, the integrated dissipation over 2000 s was observed to increase monotonically with both M and N_fil.  Quantitatively, the integrated dissipation depends on the choice of initial proportions, however we found that the shape of the dependence on N_fil and M is largely independent of initial proportions (data not shown).  Total dissipation increases with M simply because more actin is available to hydrolyze ATP during the approach to steady-state.  For large amounts of actin, increasing N_fil, the number concentration of filaments, allows increased polymerization of G^T, which constitutes a large contribution to total dissipation during the early stages of the trajectory.  Increasing N_fil also shifts the steady-state concentration of G^T downward (SI Figure <ref>), implying that more G^T has been polymerized during the approach to steady-state.  A loose analogy can be made of a one lane road compared to a multi-lane highway during heavy traffic to describe this situation.  As N_fil is increased with M fixed, the steady-state dissipation rate increases concavely, as shown in SI Figure <ref>.  The steady-state concentration of F^T also increases concavely (SI Figure <ref>), representing higher rates of ATP hydrolysis.  The contributions of each reaction to the total dissipation rate at steady-state as N_fil is varied is illustrated in SI Figure  <ref>.  The lack of dependence of the steady-state dissipation rate on M can be explained by the fact that increasing M increases the steady-state concentration of only F^D, not of any other species <cit.>.  In other words, all extra actin accumulates in the form F^D as M is increased.  This species is essentially inert, since the depolymerization rate of F^D is controlled by the concentration of filaments N_fil.  Thus the steady-state dissipation has no dependence on the total amount of actin.  Furthermore, it has no dependence on the initial concentrations, since it is known that the steady-state vector of concentrations does not depend on initial conditions <cit.>.  
	
	

		
	

	
	

		
	

	
	

		
	

	
	



§ DATA ACCESSIBILITY

MEDYAN is available under copyright for download at http://www.medyan.org/.  The version code and data analyzed here without mechanochemical feedback are available through the digital repository at the University of Maryland:  https://doi.org/10.13016/t8oa-9qra.  


§ AUTHORS' CONTRIBUTIONS

G.P. conceived the project. G.P. and C.J. consulted on all aspects of the project.  C.F. developed and implemented the new MEDYAN code, performed the simulations, analyzed results, and wrote the manuscript.  All authors contributed to the mean-field model, edited the manuscript, and approved it for submission.  


§ COMPETING INTERESTS

The authors declare no competing or financial interests.  


§ FUNDING

This work was supported by the following grants from the National Science Foundation: NSF 1632976, NSF DMR-1506969, and NSF CHE-1800418.


§ ACKNOWLEDGMENTS

The authors gratefully acknowledge A. Chandresekaran, Q. Ni, and  H. Ni for their improvements of the manuscript and helpful discussions.  



unsrt

10
	
	fletcher2010cell
	Daniel A Fletcher and R Dyche Mullins.
	Cell mechanics and the cytoskeleton.
	Nature, 463(7280):485, 2010.
	
	delon2007integrins
	Isabelle Delon and Nicholas H Brown.
	Integrins and the actin cytoskeleton.
	Current opinion in cell biology, 19(1):43–50, 2007.
	
	burlacu1992distribution
	S Burlacu, PA Janmey, and J Borejdo.
	Distribution of actin filament lengths measured by fluorescence
	microscopy.
	American Journal of Physiology-Cell Physiology,
	262(3):C569–C577, 1992.
	
	fritzsche2016actin
	Marco Fritzsche, Christoph Erlenkämper, Emad Moeendarbary, Guillaume
	Charras, and Karsten Kruse.
	Actin kinetics shapes cortical network structure and mechanics.
	Science advances, 2(4):e1501337, 2016.
	
	stricker2010mechanics
	Jonathan Stricker, Tobias Falzone, and Margaret L Gardel.
	Mechanics of the f-actin cytoskeleton.
	Journal of biomechanics, 43(1):9–14, 2010.
	
	wachsstock1993affinity
	Daniel H Wachsstock, WH Schwartz, and Thomas D Pollard.
	Affinity of alpha-actinin for actin determines the structure and
	mechanical properties of actin filament gels.
	Biophysical journal, 65(1):205–214, 1993.
	
	keller2000mechanochemistry
	David Keller and Carlos Bustamante.
	The mechanochemistry of molecular motors.
	Biophysical Journal, 78(2):541–556, 2000.
	
	albert2014stochastic
	Philipp J Albert, Thorsten Erdmann, and Ulrich S Schwarz.
	Stochastic dynamics and mechanosensitivity of myosin ii
	minifilaments.
	New Journal of Physics, 16(9):093019, 2014.
	
	alvarado2017force
	José Alvarado, Michael Sheinman, Abhinav Sharma, Fred C MacKintosh, and
	Gijsje H Koenderink.
	Force percolation of contractile active gels.
	Soft matter, 13(34):5624–5644, 2017.
	
	kovacs2003functional
	Mihály Kovács, Fei Wang, Aihua Hu, Yue Zhang, and James R Sellers.
	Functional divergence of human cytoplasmic myosin ii kinetic
	characterization of the non-muscle iia isoform.
	Journal of Biological Chemistry, 278(40):38132–38140, 2003.
	
	zhang2016thermodynamic
	Xuejun C Zhang and Wei Feng.
	Thermodynamic aspects of atp hydrolysis of actomyosin complex.
	Biophysics reports, 2(5-6):87–94, 2016.
	
	erdmann2016sensitivity
	Thorsten Erdmann, Kathrin Bartelheimer, and Ulrich S Schwarz.
	Sensitivity of small myosin ii ensembles from different isoforms to
	mechanical load and atp concentration.
	Physical Review E, 94(5):052403, 2016.
	
	mccullagh2014unraveling
	Martin McCullagh, Marissa G Saunders, and Gregory A Voth.
	Unraveling the mystery of atp hydrolysis in actin filaments.
	Journal of the American Chemical Society, 136(37):13053–13058,
	2014.
	
	li2009actin
	Xin Li, Jan Kierfeld, and Reinhard Lipowsky.
	Actin polymerization and depolymerization coupled to cooperative
	hydrolysis.
	Physical review letters, 103(4):048102, 2009.
	
	ramaswamy2010mechanics
	Sriram Ramaswamy.
	The mechanics and statistics of active matter.
	Annual Review of Condensed Matter Physics, 1:323–345, 2010.
	
	marchetti2013hydrodynamics
	M Cristina Marchetti, Jean-François Joanny, Sriram Ramaswamy,
	Tanniemola B Liverpool, Jacques Prost, Madan Rao, and R Aditi Simha.
	Hydrodynamics of soft active matter.
	Reviews of Modern Physics, 85(3):1143, 2013.
	
	hill2004free
	Terrell L Hill.
	Free energy transduction and biochemical cycle kinetics.
	Courier Corporation, 2004.
	
	tan2016self
	Tzer Han Tan, Maya Malik Garbi, Enas Abu-Shah, Junang Li, Abhinav Sharma,
	Fred C MacKintosh, Kinneret Keren, Christoph F Schmidt, and Nikta Fakhri.
	Self-organization of stress patterns drives state transitions in
	actin cortices.
	arXiv preprint arXiv:1603.07600, 2016.
	
	kohler2011structure
	Simone Köhler, Volker Schaller, and Andreas R Bausch.
	Structure formation in active networks.
	Nature materials, 10(6):462, 2011.
	
	schaller2010polar
	Volker Schaller, Christoph Weber, Christine Semmrich, Erwin Frey, and Andreas R
	Bausch.
	Polar patterns of driven filaments.
	Nature, 467(7311):73, 2010.
	
	alvarado2013molecular
	José Alvarado, Michael Sheinman, Abhinav Sharma, Fred C MacKintosh, and
	Gijsje H Koenderink.
	Molecular motors robustly drive active gels to a critically connected
	state.
	Nature Physics, 9(9):591, 2013.
	
	prigogine1967introduction
	Ilya Prigogine.
	Introduction to thermodynamics of irreversible processes.
	New York: Interscience, 1967, 3rd ed., 1967.
	
	prigogine1971biological
	Ilya Prigogine and Grégoire Nicolis.
	Biological order, structure and instabilities.
	Quarterly Reviews of Biophysics, 4(2-3):107–148, 1971.
	
	england2013statistical
	Jeremy L England.
	Statistical physics of self-replication.
	The Journal of chemical physics, 139(12):09B623_1, 2013.
	
	england2015dissipative
	Jeremy L England.
	Dissipative adaptation in driven self-assembly.
	Nature nanotechnology, 10(11):919, 2015.
	
	perunov2016statistical
	Nikolay Perunov, Robert A Marsland, and Jeremy L England.
	Statistical physics of adaptation.
	Physical Review X, 6(2):021036, 2016.
	
	kachman2017self
	Tal Kachman, Jeremy A Owen, and Jeremy L England.
	Self-organized resonance during search of a diverse chemical space.
	Physical review letters, 119(3):038001, 2017.
	
	horowitz2017spontaneous
	Jordan M Horowitz and Jeremy L England.
	Spontaneous fine-tuning to environment in many-species chemical
	reaction networks.
	Proceedings of the National Academy of Sciences,
	114(29):7565–7570, 2017.
	
	baiesi2018life
	Marco Baiesi and Christian Maes.
	Life efficiency does not always increase with the dissipation rate.
	Journal of Physics Communications, 2(4):045017, 2018.
	
	balland2005dissipative
	Martial Balland, Alain Richert, and François Gallet.
	The dissipative contribution of myosin ii in the cytoskeleton
	dynamics of myoblasts.
	European Biophysics Journal, 34(3):255–261, 2005.
	
	mitrossilis2009single
	Demosthene Mitrossilis, Jonathan Fouchard, Axel Guiroy, Nicolas Desprat,
	Nicolas Rodriguez, Ben Fabry, and Atef Asnacios.
	Single-cell response to stiffness exhibits muscle-like behavior.
	Proceedings of the National Academy of Sciences,
	106(43):18243–18248, 2009.
	
	hoffman2006consensus
	Brenton D Hoffman, Gladys Massiera, Kathleen M Van Citters, and John C Crocker.
	The consensus mechanics of cultured mammalian cells.
	Proceedings of the National Academy of Sciences,
	103(27):10259–10264, 2006.
	
	kurzawa2017dissipation
	Laetitia Kurzawa, Benoit Vianay, Fabrice Senger, Timothée Vignaud, Laurent
	Blanchoin, and Manuel Théry.
	Dissipation of contractile forces: the missing piece in cell
	mechanics.
	Molecular biology of the cell, 28(14):1825–1832, 2017.
	
	popov2016medyan
	Konstantin Popov, James Komianos, and Garegin A Papoian.
	Medyan: mechanochemical simulations of contraction and polarity
	alignment in actomyosin networks.
	PLoS computational biology, 12(4):e1004877, 2016.
	
	hu2010mechano
	Longhua Hu and Garegin A Papoian.
	Mechano-chemical feedbacks regulate actin mesh growth in
	lamellipodial protrusions.
	Biophysical journal, 98(8):1375–1384, 2010.
	
	komianos2018stochastic
	James E Komianos and Garegin A Papoian.
	Stochastic ratcheting on a funneled energy landscape is necessary for
	highly efficient contractility of actomyosin force dipoles.
	Physical Review X, 8(2):021006, 2018.
	
	ragazzon2018energy
	Giulio Ragazzon and Leonard J Prins.
	Energy consumption in chemical fuel-driven self-assembly.
	Nature nanotechnology, page 1, 2018.
	
	2019arXiv190110520F
	Carlos Floyd, Garegin A. Papoian, and Christopher Jarzynski.
	A Discrete Approximation to Gibbs Free Energy of Chemical Reactions
		is Needed for Accurately Calculating Entropy Production in Mesoscopic
		Simulations.
	arXiv e-prints, page arXiv:1901.10520, January 2019.
	
	brooks2009nonequilibrium
	FJ Brooks and AE Carlsson.
	Nonequilibrium actin polymerization treated by a truncated
	rate-equation method.
	Physical Review E, 79(3):031914, 2009.
	
	floyd2017low
	Carlos Floyd, Christopher Jarzynski, and Garegin Papoian.
	Low-dimensional manifold of actin polymerization dynamics.
	New Journal of Physics, 19(12):125012, 2017.
	
	leys2013detecting
	Christophe Leys, Christophe Ley, Olivier Klein, Philippe Bernard, and Laurent
	Licata.
	Detecting outliers: Do not use standard deviation around the mean,
	use absolute deviation around the median.
	Journal of Experimental Social Psychology, 49(4):764–766,
	2013.
	
	huber2011robust
	Peter J Huber.
	Robust statistics.
	Springer, 2011.
	
	erban2014multiscale
	Radek Erban, Mark B Flegg, and Garegin A Papoian.
	Multiscale stochastic reaction–diffusion modeling: application to
	actin dynamics in filopodia.
	Bulletin of mathematical biology, 76(4):799–818, 2014.
	
	hu2013molecular
	Longhua Hu and Garegin A Papoian.
	Molecular transport modulates the adaptive response of branched actin
	networks to an external force.
	The Journal of Physical Chemistry B, 117(42):13388–13396,
	2013.
	
	jones2014scipy
	Eric Jones, Travis Oliphant, and Pearu Peterson.
	{SciPy}: open source scientific tools for {Python}.
	2014.
	
	limpert2001log
	Eckhard Limpert, Werner A Stahel, and Markus Abbt.
	Log-normal distributions across the sciences: Keys and clues: On the
	charms of statistics, and how mechanical models resembling gambling machines
	offer a link to a handy way to characterize log-normal distributions, which
	can provide deeper insight into variability and probability—normal or
	log-normal: That is the question.
	AIBS Bulletin, 51(5):341–352, 2001.
	
	furusawa2005ubiquity
	Chikara Furusawa, Takao Suzuki, Akiko Kashiwagi, Tetsuya Yomo, and Kunihiko
	Kaneko.
	Ubiquity of log-normal distributions in intra-cellular reaction
	dynamics.
	Biophysics, 1:25–31, 2005.
	
	kundu2005discriminating
	Debasis Kundu and Anubhav Manglick.
	Discriminating between the log-normal and gamma distributions.
	Journal of the Applied Statistical Sciences, 14:175–187, 2005.
	
	bak1988self
	Per Bak, Chao Tang, and Kurt Wiesenfeld.
	Self-organized criticality.
	Physical review A, 38(1):364, 1988.
	
	bak1989earthquakes
	Per Bak and Chao Tang.
	Earthquakes as a self-organized critical phenomenon.
	Journal of Geophysical Research: Solid Earth,
	94(B11):15635–15637, 1989.
	
	sornette1989self
	A Sornette and D Sornette.
	Self-organized criticality and earthquakes.
	EPL (Europhysics Letters), 9(3):197, 1989.
	
	bottiglieri2007off
	M Bottiglieri and C Godano.
	On-off intermittency in earthquake occurrence.
	Physical Review E, 75(2):026101, 2007.
	
	cardamone2011cytoskeletal
	Luca Cardamone, Alessandro Laio, Vincent Torre, Rajesh Shahapure, and Antonio
	DeSimone.
	Cytoskeletal actin networks in motile cells are critically
	self-organized systems synchronized by mechanical interactions.
	Proceedings of the National Academy of Sciences,
	108(34):13978–13983, 2011.
	
	freedman2018nonequilibrium
	Simon L Freedman, Glen M Hocky, Shiladitya Banerjee, and Aaron R Dinner.
	Nonequilibrium phase diagrams for actomyosin networks.
	Soft matter, 14(37):7740–7747, 2018.
	
	seara2018entropy
	Daniel S Seara, Vikrant Yadav, Ian Linsmeier, A Pasha Tabatabai, Patrick W
	Oakes, SM Ali Tabei, Shiladitya Banerjee, and Michael P Murrell.
	Entropy production rate is maximized in non-contractile actomyosin.
	Nature communications, 9(1):4948, 2018.
	
	battle2016broken
	Christopher Battle, Chase P Broedersz, Nikta Fakhri, Veikko F Geyer, Jonathon
	Howard, Christoph F Schmidt, and Fred C MacKintosh.
	Broken detailed balance at mesoscopic scales in active biological
	systems.
	Science, 352(6285):604–607, 2016.
	
	falzone2015entangled
	Tobias T Falzone, Savanna Blair, and Rae M Robertson-Anderson.
	Entangled f-actin displays a unique crossover to microscale
	nonlinearity dominated by entanglement segment dynamics.
	Soft matter, 11(22):4418–4423, 2015.
	
	fujiwara2007polymerization
	Ikuko Fujiwara, Dimitrios Vavylonis, and Thomas D Pollard.
	Polymerization kinetics of adp-and adp-pi-actin determined by
	fluorescence microscopy.
	Proceedings of the National Academy of Sciences,
	104(21):8827–8832, 2007.
	
	mcfadden2017filament
	William M McFadden, Patrick M McCall, Margaret L Gardel, and Edwin M Munro.
	Filament turnover tunes both force generation and dissipation to
	control long-range flows in a model actomyosin cortex.
	PLoS computational biology, 13(12):e1005811, 2017.
	
	duda2009thermodynamics
	Andrzej Duda and Adam Kowalski.
	Thermodynamics and kinetics of ring-opening polymerization.
	Handbook of ring-opening polymerization, pages 1–51, 2009.
	
	erdmann2013stochastic
	Thorsten Erdmann, Philipp J Albert, and Ulrich S Schwarz.
	Stochastic dynamics of small ensembles of non-processive molecular
	motors: The parallel cluster model.
	The Journal of chemical physics, 139(17):11B604_1, 2013.
	
	howard2001mechanics
	Jonathon Howard.
	Mechanics of motor proteins and the cytoskeleton.
	Sinauer associates Sunderland, MA, 2001.
	
	sakamoto2008direct
	Takeshi Sakamoto, Martin R Webb, Eva Forgacs, Howard D White, and James R
	Sellers.
	Direct observation of the mechanochemical coupling in myosin va
	during processive movement.
	Nature, 455(7209):128, 2008.
	
	dufort1996profilin
	Paul A Dufort and Charles J Lumsden.
	How profilin/barbed-end synergy controls actin polymerization: A
	kinetic model of the atp hydrolysis circuit.
	Cytoskeleton, 35(4):309–330, 1996.
	
	yarmola2008effect
	Elena G Yarmola, Dmitri A Dranishnikov, and Michael R Bubb.
	Effect of profilin on actin critical concentration: a theoretical
	analysis.
	Biophysical journal, 95(12):5544–5573, 2008.
	
	nowak1988kinetics
	Ewa Nowak and Roger S Goody.
	Kinetics of adenosine 5'-triphosphate and adenosine 5'-diphosphate
	interaction with g-actin.
	Biochemistry, 27(23):8613–8617, 1988.
	
	kinosian1993nucleotide
	HJ Kinosian, LA Selden, JE Estes, and LC Gershman.
	Nucleotide binding to actin. cation dependence of nucleotide
	dissociation and exchange rates.
	Journal of Biological Chemistry, 268(12):8683–8691, 1993.
	
	selden1999impact
	Lynn A Selden, Henry J Kinosian, James E Estes, and Lewis C Gershman.
	Impact of profilin on actin-bound nucleotide exchange and actin
	polymerization dynamics.
	Biochemistry, 38(9):2769–2778, 1999.
	
	milo2015cell
	Ron Milo and Rob Phillips.
	Cell biology by the numbers.
	Garland Science, 2015.
	
	nicolis1977
	G Nicolis and I Prigogine.
	I (1977) self-organization in nonequilibrium systems.
	From Dissipative structures to Order through Fluctuations, Mir,
		Moscow, Russia.
	
	kondepudi2014modern
	Dilip Kondepudi and Ilya Prigogine.
	Modern thermodynamics: from heat engines to dissipative
		structures.
	John Wiley & Sons, 2014.
	
	glansdorff1954proprietes
	P Glansdorff and I Prigogine.
	Sur les propriétés différentielles de la production
	d'entropie.
	Physica, 20(7-12):773–780, 1954.
	
	peskin1993cellular
	Charles S Peskin, Garrett M Odell, and George F Oster.
	Cellular motions and thermal fluctuations: the brownian ratchet.
	Biophysical journal, 65(1):316–324, 1993.
	








