Insights from Graph Theory on the Morphologies of Actomyosin Networks with Multilinkers
Yossi Eliaz, Francois Nedelec, Greg Morrison, Herbert Levine, Margaret S. Cheung
11 Multivalent actin-binding proteins augment the variety of morphologies in actomyosin networks
Yossi Eliaz , Francois Nedelec , Greg Morrison , Herbert Levine , and Margaret S. Cheung * Department of Physics, University of Houston, Houston, Texas 77204, USA Center for Theoretical Biological Physics, Rice University, Houston, Texas 77005, USA Sainsbury Laboratory, Cambridge University, Bateman Street, CB2 1LR Cambridge, UK Department of Physics, Northeastern University, Boston, Massachusetts 02115, USA Department of Bioengineering, Rice University, Houston, Texas 77030, USA *Corresponding author: [email protected]
ABSTRACT
We study the role of multivalent actin-binding proteins in remodeling actomyosin networks, using a numerical model of semiflexible filaments. We find that the valency of actin-binding proteins greatly affects the organization of actomyosin networks. The rich variety of emergent architectures and morphologies inspired us to develop order parameters based on network theory, to quantitatively characterize the self-assembled organization in the simulated network. First, actin-binding proteins with a valency greater than two promote filament bundles and large filament clusters to a much greater extent than bivalent multilinkers. Second, active myosin-like motor proteins promote rich phases of actomyosin, as characterized by network theory. Lastly, we find that both high content of multilinkers and active motors are necessary to the reduction of node degree correlation in the network of filaments. Our work provides a framework to describe the dendritic structures of actomyosin networks mediated by ABP with the ability of binding multiple filaments, such as Calcium/Calmodulin-dependent Kinase II (CaMKII). It establishes a theoretical basis to promote the use of network theory in analyzing experiments, and help us to understand the role of dendritic spine growth in synaptic plasticity.
I. INTRODUCTION
Living cells actively regulate the morphologies of actomyosin networks to control the force produced by these networks during various cellular processes [1, 2]. They do so, in part, by controlling the activity of specific actin-binding proteins (ABPs) that either bind to actin monomers, the filaments, or both. Some of the ABPs, such as Myosin II (a molecular motor) [3] or CaMKII [4, 5], require the hydrolysis of ATP to actively transmit mechanical changes in actomyosin networks, whereas passive ABPs such as α -actinin coordinate contractions and control the density of actomyosin networks by cross-linking thin filaments [6]. Motors and crosslinkers collectively regulate the morphology of a dendritic spine in neuronal cells as well as the movement of a cell [7, 8] by altering the connectivity of the actin filaments, promoting the formation of a biologically functional architecture [4, 9]. There have been extensive studies on ABPs in maintaining or modifying actomyosin structures. However, most ABPs and actin-related proteins (ARPs), such as Arp2/3 [10, 11], operate on two filaments as “bivalent linkers”. ABPs with valency greater than 2 exist. The Tau protein is a multivalent linker [12] and neurodegeneration mediator [13]. The Ca /calmodulin-dependent protein kinase II (CaMKII) protein is also a multivalent actin-binding protein [4, 5], which plays a mechanobiological role in tuning the intricate shape of dendritic spines in neurons. We hypothesize that with their ability to bind more than two filaments, these multi-linking ABPs, or “multilinkers”, provide a tunable feature to diversify the morphology of actomyosin networks. We tested our hypothesis by using Cytosim [14], a software developed to simulate mesoscopic cytoskeletal biological systems [9, 15]. In this system, polymeric actin is represented as filaments with bending elasticity and linkers as stochastic entities that can bind and connect these filaments into a network. Extending Cytosim’s modeling of motors, actin filaments and bivalent actin-binding proteins, we incorporate a model for multivalent ABPs. With this addition, Cytosim allowed us to investigate the structural and dynamical principles that govern the organization of a random actomyosin networks, depending on motor concentration and on the valency of the multivalent-ABPs. We have discovered that the multivalent ABPs enriched the variety in the dynamics and structures of actomyosin networks as they evolve heterogeneously in space with time. Intrinsic order parameters for describing the phase diagram of a single- component matter are insufficient to fully grasp the comprehensive states of dendritic actomyosins in a wide parameter space. There is a need for reduced representation of the networks in terms of simple order parameters, yet still sophisticated enough to capture the hierarchical features in the actomyosin networks mediated by distinct ABPs. Some of the ABPs such as CaMKII are directly coupled to calcium signaling pathways [4], which provides a unique timed control of actomyosin reorganization activated by external stimuli through oscillatory variations in calcium signals. We compared the order parameters from different domains of science in terms of polymer physics, gelation, and graph theory to characterize and interpret the simulation results. We ranked these order parameters by investigating their correlations. We found that the order parameters from network theory capture local and non-local features in a dendritic actomyosin network. We established the state diagrams in plane of selected network theory order parameters to identify and quantify the effects of motor activity and multilinker concentration on actomyosin architecture. We have discovered subtle phases in dendritic growths of actomyosin in some parameter range of a state diagram that permits experimental validations.
II. MOTIVATION, BIOPHYSICAL BACKGROUND, AND HYPOTHESIS A. Dendritic spines in neurons
Neurons are brain cells that communicate using electrical and chemical signals through action potential signaling and neurotransmitter release. Two neurons are “in touch” without physically touching one another. At the junction between the two neurons, called a synapse, the presynaptic neuron signals the postsynaptic one. Dendritic spines are protrusions rising from postsynaptic neuron which form the receiver end of the synapse. Spine volume ranges between −3 μ m and μ m [16], and each neuron may grow thousands of spines. Consistent with the Hebbian theory of learning and memory, spines expand during long-term potentiation (LTP) and shrink during long-term depression (LTD) [17]. These learning processes do not solely involve chemical and electrical signals but also cellular mechanics. Experiments of high-resolution imaging have shown that LTP induction causes morphological changes in dendritic spines [18, 19]. Besides the shrinking or expanding during LTD and LTP, the generation and destruction of spines have been observed [16]. The spine is rich with ABPs such as CaMKII crucial for maintaining the information and structural of a spine [20, 21]. B. CaMKII, morphogenesis, and synaptic plasticity
The Ca /calmodulin-dependent kinase II [22] (CaMKII) holoenzyme protein serves a dual role as a Ca signaling decoder and as a structural agent in directing calcium signals to change the makeup of actomyosin networks in a spine. It accounts for up to 2% of the mammalian brain proteome [22, 23]. In vertebrates, there are four major CaMKII isoforms observed in over 40 different splice variants. Overall, these variants are genetically expressed in diverse tissue types. The two most prevalent isoforms in the brain are CaMKII α and CaMKII β [5, 24] and the number of the protein domains in these multimeric complexes varies from 12 to 14. Each unit is made up of several connected domains that serve distinctive functions [4]. Besides being a biochemical catalyzer and having autophosphorylation capabilities [25] in one domains, each monomeric unit is capable of binding actin filaments. To date, the highest actin-binding valency observed for CaMKII is six for the dodecamer (12-mer) CaMKII β isoform [5]. This makes CaMKII a rare kind of multivalent actin-binding protein [4, 5, 26]. It motivates the development of our hypothesis that ABP’s valency affects the general morphology of actomyosin networks. It may elucidate a functional role for the multivalent association nature of the CaMKII, which regulates postsynaptic computational resources by changing the morphology of the dendritic spine [27, 28]. CaMKII is essential for long-term potentiation [29] in a persistent strengthening of synapses. The morphological changes over time is called synaptic plasticity [4]. C. Harnessing network theory to study structural biophysical systems
Contact maps are 2D representation of a 3D fold of a linear object (e.g. a polymer). With this approach, physical proximity is represented by an element of a 2D matrix where both columns and lines represent the linear position along the object (e.g. the residue number). Contact maps are widely used in structural biology, from simulations of protein folding [30, 31] to the visualization of the human genome’s at one thousand base-pair resolution [32]. Nonetheless, it remains difficult to exploit such maps and more generally the formalism of network and graph theory to investigate mesoscopic biological systems [33]. Although it is possible to label individual actin filaments or motor proteins, the resolution of light microscopy limits the ability to distinguish between individual filaments and decipher their linkages. The advent of superresolution 3D imaging techniques may enable direct observation of individual filaments and the connections between them in vivo , permitting a more complete picture of the global topology of these systems [34] [35]. It is thus essential to develop the tools to better understand not only the microscopic details underlying the filaments, crosslinks, or motors using traditional biophysical modeling, but also to understand how the mesoscopic network structure, through non-local order parameters, relates to biological function. The applications of network approaches are broad in physics: the physics of large networks [33, 34], renormalization group of tensor networks [36], quantum entanglement [37], and the spread of epidemic diseases [38, 39]. In this paper we will apply it in the actomyosin networks by bridge a gap between detailed dynamics and structure of actomyosin networks from simulations and observables in high-resolution experiments. We will survey and compare order parameters widely employed from polymer physics and soft matter, as well as the order parameters arising from graph theory on the growth of dendritic actomyosin networks from Cytosim. Our far reaching goal is to leverage graph theory to develop a statistical physics framework of nonequilibrium ensembles and landscape reconstruction [40, 41] for actomyosin networks.
III. MODELS AND METHODS A. Coarse-grained model of actomyosin networks
We used a coarse-grained model to study the morphology and structure of actomyosin networks. In our 3D model, actin filaments are represented as incompressible bendable polar fibers with rigidity of 0.075 pN μ m [42]. Myosin motors have two walking heads, each of which operates independently and walks toward a filament’s plus end. Additionally, we have two types of linker species: one that resembles an α -actinin bivalent crosslinker and second, the multivalent linker inspired by CaMKII with varied valencies referred to as a 'multilinker' in this article. As depicted in Fig. 1, the simulation contains filaments, passive bivalent crosslinkers, passive multivalent linkers (multilinkers), and active bivalent motors. B. Collective Constrained Langevin dynamics
The vector 𝐱𝐱 (t) contains the coordinates of the N three-dimensional vertices describing the physical objects in the system at time t . For a fixed temperature 𝑇𝑇 , the equation of motion is written in a Langevin form: 𝑑𝑑𝐱𝐱 (t) = 𝛍𝛍 𝐟𝐟 𝐭𝐭𝐭𝐭𝐭𝐭 ( 𝐱𝐱 , t) 𝑑𝑑 t + 𝑑𝑑𝐁𝐁 (t), where 𝛍𝛍 is a
3N × 3N diagonal matrix consisting of the mobility coefficients of all objects, 𝐟𝐟 𝐭𝐭𝐭𝐭𝐭𝐭 ( 𝐱𝐱 , t) is a vector contains all the forces acting on the points 𝐱𝐱 (t) at time t. 𝐁𝐁 (t) recapitulates the random Brownian noise vector due to molecular collisions and its 𝑖𝑖 th coordinate, 𝐁𝐁 i (t) which is a temperature-dependent value drawn from a normal distribution with a zero mean and a standard deviation equals to � 𝑑𝑑𝑑𝑑𝐃𝐃 i , where 𝑑𝑑𝑑𝑑 is the integration time. Per Einstein’s relation we set the diffusion coefficient of the 𝑖𝑖 th random molecular degree of freedom to be 𝐃𝐃 i ≡ 𝛍𝛍 ii 𝑘𝑘 B 𝑇𝑇 , where 𝑘𝑘 B is the Boltzmann constant. B. Motor activity and binding and unbinding events
We used the software Cytosim [14] to propagate the equation of motions. After the collective Brownian mechanics has been calculated, Cytosim executes two sub-routines to account for chemical processes such as binding, unbinding, and motor walking. The first subroutine simulates binding and unbinding events of actin-binding proteins according to the 𝑘𝑘 on and 𝑘𝑘 off rates. The second subroutine emulates motor activity, computes the motor’s exerted forces on the filaments and recalculates the locations of motors on the filaments [14]. C. Shape and mobility of fibers and solid objects
Filaments and multilinkers are described by vertices with three spatial coordinates each. A filament is an elastic, polar, and non-extendible rod represented by a set of 𝑝𝑝 equidistant vertices { 𝑚𝑚 𝑖𝑖 } 𝑖𝑖=1𝑝𝑝 , where 𝑚𝑚 and 𝑚𝑚 𝑝𝑝 are the minus and the plus ends, respectively. A solid object in Cytosim is a nondeformable set of vertices with fixed size and shape. Filaments emulate filamentous actin, whereas the skeleton of the multilinkers is a solid object in Cytosim. A solid object is a set of 𝑘𝑘 vertices { 𝑠𝑠 𝑖𝑖 } 𝑖𝑖=1𝑘𝑘 with hard constraints on their pairwise distances ( | 𝑠𝑠 𝑖𝑖 − 𝑠𝑠 𝑗𝑗 | = 𝑑𝑑 𝑖𝑖𝑗𝑗 = 𝑑𝑑 𝑗𝑗𝑖𝑖 ); e.g., the multilinker’s binding entities that can possibly bind to filaments lay on the surface of a solid sphere, which effectively only has translational and rotational degrees of freedom. D. Interactions between objects
In Cytosim, any interaction between objects is expressed by a linearization 𝐟𝐟 ( 𝐫𝐫 , t) = 𝐀𝐀 (t) 𝐫𝐫 + 𝐠𝐠 (t) . The matrix 𝐀𝐀 (t) and the vector 𝐠𝐠 (t) contain the contributions from all the elementary interactions limited to the constant and linear terms of the Taylor expansion. Two points 𝐫𝐫 𝑎𝑎 and 𝐫𝐫 𝑏𝑏 from different objects ( 𝑎𝑎 and 𝑏𝑏 ) can be connected by a link with Hookean stiffness 𝑘𝑘 . The forces between the points read: 𝐟𝐟 𝑎𝑎 = −𝐟𝐟 𝑏𝑏 = 𝑘𝑘 � − 𝑟𝑟 | 𝐫𝐫 𝑏𝑏 – 𝐫𝐫 𝑎𝑎 | � ( 𝐫𝐫 𝑏𝑏 – 𝐫𝐫 𝑎𝑎 ) , (1) where 𝑟𝑟 ≥ is the resting length of the link. When a motor is attached to a fiber, the motor’s position is given by a curvilinear abscissa measured from a fixed reference point on the fiber ( 𝐱𝐱 ). The abscissa is increased by 𝛿𝛿 = 𝜏𝜏𝑣𝑣 maxmotor � − 𝑓𝑓𝑓𝑓 𝑠𝑠𝑠𝑠𝑎𝑎𝑠𝑠𝑠𝑠 � , where 𝑣𝑣 maxmotor is a constant real value representing the maximal speed of a motor walking along a fiber. Its sign dictates the motor’s tendency to walk toward filament’s minus or plus end. The value of 𝑓𝑓 is the load that the motor experiences projected along the direction of the filament on which the motor walks. The stall force 𝑓𝑓 𝑠𝑠𝑠𝑠𝑎𝑎𝑠𝑠𝑠𝑠 is the amount of force that is sufficient to stop the motor from moving. In general, Cytosim uses a force-dependent unbinding rate 𝑘𝑘 off = 𝑘𝑘 exp � | 𝑓𝑓 | 𝑓𝑓 � to model dissociation from the fiber. Motors have constant binding 𝑘𝑘 onmotor and unbinding rates 𝑘𝑘 offmotor . Namely, they work in the limit of constant unbinding rate ( 𝑓𝑓 → ∞ ). In addition, we set 𝑣𝑣 maxmotor = 0.2 μ m s −1 and 𝑓𝑓 𝑠𝑠𝑠𝑠𝑎𝑎𝑠𝑠𝑠𝑠 = 6 pN . The bivalent crosslinker and multilinkers have the unbinding and binding rates: 𝑘𝑘 onlinker and 𝑘𝑘 offlinker . The interaction of a linker and filaments is assumed to be Hookean with zero resting length and a stiffness of 50 pN μ m −1 . FIG. 1. Illustration of Cytosim’s model and its components. A vertex-based model of actin filament in (a) with an accompanying crystal structure of 10 globular actin (G-actin) monomers build up a filamentous actin (F-actin) in (d), based on PDB ID: 3J8I. In (b) and (c) are the side and front views of CaMKII crystal structure (PDB ID: 3SOA) visualizing only the association domains of CaMKII. The protein complexes in (b), (c), and (d) were generated using UCSF Chimera [43]. A multilinker of 𝑣𝑣 -valency was built using a sphere with 𝑣𝑣 active sites arranged on its surface such as to maximize the minimum distance between any two binding entities. The special case of tetravalent ( 𝑣𝑣 = 4 ) linker with a radius attached to two filaments is depicted here (e). An α -actinin-like bivalent crosslinker is modeled as a spring between two filaments as shown in (f). In (g) and (h), there is an illustration of binding and unbinding of a diffusive particle such as: a motor, a crosslinker, or a multilinker from a filament. A myosin-like motor, walking on a filament and exerting a force, is depicted in (i) with the motor’s force-velocity relationship. CaMKII’s functional domains are shown in (j) including the association domain that binds F-actin [4] and the catalytic domain that binds to ATP and is regulated by the regulatory (Reg.) domain, which becomes active by calcium/calmodulin signaling [44]. E. Excluded Volume Interactions Between Multilinkers
Cytosim incorporates excluded volume and “steric” effects between objects. Generally, it supports both attractive and repulsive forces as the interaction is taken as piecewise linear force: 𝑓𝑓 ( 𝑑𝑑 ) = �𝑘𝑘 + ( 𝑑𝑑 − 𝑑𝑑 ), 𝑑𝑑 ≤ 𝑑𝑑 𝑘𝑘 − ( 𝑑𝑑 − 𝑑𝑑 ), 𝑑𝑑 < 𝑑𝑑 < 𝑑𝑑 + 𝑑𝑑 𝑑𝑑 + 𝑑𝑑 ≤ 𝑑𝑑 , (2) where 𝑑𝑑 is the distance between two interacting elements, 𝑑𝑑 is the distance at which the elements are at equilibrium, and 𝑑𝑑 is an additional distance that confers an attractive component to the interaction in the range 𝑑𝑑 < 𝑑𝑑 < 𝑑𝑑 + 𝑑𝑑 . We use only a repulsive steric effect between multilinkers with 𝑑𝑑 = 5 nm , 𝑘𝑘 − = 0 , and 𝑘𝑘 + = 500 pN μ m −1 and thus 𝑑𝑑 = 0 . F. General Simulation’s Setup Parameters
Our system is a m cube with a viscosity of μ m −2 s . The filaments have length of μ m and are uniformly distributed in the cubic system. The actin filaments have a flexural rigidity of μ m (corresponding to a persistence length ~ m at the configuration working temperature, 𝑘𝑘 B 𝑇𝑇 = 4.2 pNnm ) [42]. An α -actinin-like bivalent crosslinker is modeled as a Hookean spring of zero resting length between two binding entities with spring stiffness 𝑘𝑘 α−actinin = 250 pN μ m −1 . Each binding entity has a binding range of , binding rate 𝑘𝑘 onlinker = 5 s −1 , and unbinding rate 𝑘𝑘 offlinker = 0.1 s −1 . A myosin-like motor is an inextensible object with two independent motor heads that can walk on a filament, separated by a resting length of
100 nm . When attached, the motor and filament are Hookean with a stiffness of 𝑘𝑘 motor =250 pN μ m −1 . Motors unbind at the same rate as linkers, 𝑘𝑘 offmotor = 0.1 s −1 , and bind twice as fast as linkers, 𝑘𝑘 onmotor = 10 s −1 , and each motor site has a binding range of
50 nm . A multilinker of valency 𝑣𝑣 has 𝑣𝑣 hands residing on the surface of a sphere of radius and are linked to the sphere with a stiffness of 𝑘𝑘 multilinker = 200 pN μ m −1 . A multilinker’s binding entity has a binding range of and binding and unbinding rates equal to the rates of a crosslinker. The parameters systematically varied in the simulations are the total number of multilinkers in the system ( 𝑁𝑁 multilinkers ∈ {250,500,1000} ), their valencies ( 𝜈𝜈 ∈ {2,3,4,5,6,7} ), the number of motors ( 𝑁𝑁 motors ∈ {250,500,1000} ), the number of filaments ( 𝑁𝑁 filaments ∈ {250,500,1000} ), and the number of α -actin-like crosslinkers 𝑁𝑁 crosslinkers ∈ {10,250,500,1000} ). We explored the Cartesian product of all possible values of the above five parameters. For example, one possible system configuration included 250 pentavalent ( 𝜈𝜈 = 5) multilinkers, 500 filaments, 1000 motors, and 10 crosslinkers. For readability, Table I summarizes the possible values of parameters and the physical properties of objects in the simulation: TABLE I. System configuration parameter values for the computational model, with, when applicable, the set values that has been explored. In all simulations, these parameters are held constant from start to end, and they are only varied between simulations Parameter Description Values 𝑣𝑣 Valency of a multilinker {2,3,4,5,6,7} 𝑘𝑘 multilinker Hookean stiffness of a multilinker binding entity
200 pN μ m −1 𝑘𝑘 crosslinker Hookean stiffness of a crosslinker
250 pN μ m −1 𝑘𝑘 offlinker Unbinding rate of a multilinker or a crosslinker binding entity −1 𝑘𝑘 onlinker Binding rate of a multilinker or a crosslinker binding entity −1 𝑘𝑘 offmotor Unbinding rate of a motor binding entity −1 𝑘𝑘 onmotor Binding rate of a motor binding entity
10 s −1 𝑁𝑁 motors Number of motors {10,250,500,1000} 𝑁𝑁 filaments Number of filaments {250,500,1000} 𝑁𝑁 crosslinker Number of α -actinin crosslinkers {10,250,500,1000} V 𝑏𝑏𝑏𝑏𝑏𝑏 Volume of the simulation box m 𝑁𝑁 multilinkers Number of multilinkers {250,500,1000}
G. Order Parameters
We explore here the use of various order parameters to analyze actomyosin network simulations. An order parameter should reflect an interesting feature of the system, ideally reducing the configuration of a system characterized by a myriad of degrees of freedom into a single mesoscopic physically meaningful value. First, we examine geometrical quantities such as the radius of gyration (R g ), the shape (S), and the asphericity parameters ( ∆ ) [45] [46] among all the filament vertices. We, next use gel-sol theory order parameters, which consider the molecular connectivity of actin-binding proteins (ABP), i.e., motors and linkers, to filaments. These gel-sol order parameters measure the number of clusters, the number of filaments in any cluster, the largest cluster size, the smallest cluster size, the average cluster size, and the standard deviation in cluster sizes. Finally, we calculate network theory order parameters using undirected graph representation of the actomyosin network: 𝐆𝐆 = (V, E) . The graph has a set of vertices V , which correspond to filaments, and a set of edges E , which accounts for the connections between vertices. An edge exists between node 𝑖𝑖 and node 𝑗𝑗 if the minimal Euclidean distance between the center of mass of two segments is less than 𝑑𝑑 cutoff = 200 nm . We define the distance 𝑑𝑑 between two filaments to be the minimal distance among all possible pairwise distances of their segments’ center of masses. Thus, the graph 𝐆𝐆 is constructed only from the information of distances between filaments and without any consideration or knowledge of the microscopic connectivity between ABPs and filaments. The network-theory-based list of order parameters includes the number of communities [47], average clustering [48], clique number of the graph (the maximal clique size), average closeness centrality [49, 50], average eigenvector centrality [51], mean average neighbor degree [52], degree assortativity coefficient [53], and graph density. The Python package, NetworkX 2.2 [54], was used to construct the graph 𝐆𝐆 and to compute the order parameters.
1. Distribution of actin vertices (monomers) in space
We explored the dynamical changes in shape and structure of networks in the presence of multilinkers for both low and high motor concentrations. We used three order parameters from the protein folding field [45]: the radius of gyration 𝑅𝑅 𝑔𝑔 (the variance of all the positions of filament vertices), which is a proxy for the macroscopic structure, the shape parameter − ≤ S ≤ , and the asphericity measure ≤ Δ ≤ . The shape order parameter [45] specifies how prolate ( S > 0 ) or oblate (
S < 0 ) is the conformation of an actomyosin network, while asphericity measures how an actomyosin network conformation differs from a perfect sphere ( Δ = 0 ). The fluctuations in these order parameters reveal the stress relaxation the system undergoes through, and its response at different ABP-filaments concentration ratios. Using spectral and fluctuation analysis, we evaluate the differences [55] at high versus low motor activity. In order to define 𝑅𝑅 𝑔𝑔 , Δ , and S , we first define the moment of inertia tensor of filament vertices at time 𝑑𝑑 : 𝐓𝐓 𝛼𝛼𝛼𝛼 ( 𝑑𝑑 ) = ∑ �𝐫𝐫 𝑖𝑖𝛼𝛼 ( 𝑑𝑑 ) − 𝐫𝐫 𝑗𝑗𝛼𝛼 ( 𝑑𝑑 ) � 𝑁𝑁𝑖𝑖 , 𝑗𝑗=1 �𝐫𝐫 𝑖𝑖𝛼𝛼 ( 𝑑𝑑 ) − 𝐫𝐫 𝑗𝑗𝛼𝛼 ( 𝑑𝑑 ) � , ( ) where 𝐫𝐫 𝑖𝑖𝛼𝛼 ( 𝑑𝑑 ) is the 𝛼𝛼 -component of a filament vertex, 𝑁𝑁 is the number of filament vertices, and 𝛼𝛼 , 𝛽𝛽 ∈ {x, y, z} are the indices of the Cartesian elements. Then the radius of gyration is given by the sum of the eigenvalues, 𝜆𝜆 ≥ 𝜆𝜆 ≥ 𝜆𝜆 of 𝐓𝐓 at time 𝑑𝑑 : 𝑅𝑅 𝑔𝑔 ( 𝑑𝑑 ) = �𝑑𝑑𝑟𝑟𝐓𝐓 ( 𝑑𝑑 ) = �∑ 𝜆𝜆 𝑖𝑖 ( 𝑑𝑑 ) . ( ) We let 𝜆𝜆̅ ( 𝑑𝑑 ) = 𝑠𝑠𝑟𝑟𝐓𝐓 ( 𝑠𝑠 ) be the average tensor trace of 𝐓𝐓 ( 𝑑𝑑 ) and define the asphericity: Δ (t) =
32 �∑ �𝜆𝜆 𝑖𝑖 ( 𝑠𝑠 ) −𝜆𝜆� ( 𝑠𝑠 ) � ��𝑠𝑠𝑟𝑟𝐓𝐓 ( 𝑠𝑠 ) � . ( ) Finally, the shape parameter reads: S( 𝑑𝑑 ) = ∏ ( 𝜆𝜆 𝑖𝑖 ( 𝑠𝑠 ) −𝜆𝜆� ( 𝑠𝑠 )) � 𝑠𝑠𝑟𝑟𝐓𝐓 ( 𝑠𝑠 ) � . ( ) To calculate the power spectral density (PSD) of each order parameter, we use Welch’s method for spectral density approximation [56] with a flat top window size of 10 samples and a window overlap of . Then we compare the PSD of 𝑅𝑅 𝑔𝑔 ( 𝑑𝑑 ) , Δ ( 𝑑𝑑 ) , and S( 𝑑𝑑 ) signals in a low motor concentration versus a high motor concentration. Then, the ratio between any two power measurements, 𝑃𝑃 and 𝑃𝑃 , are computed in decibels following the convention � 𝑃𝑃 𝑃𝑃 � . The value of Δ (t) indicates the anisotropy of the system and S( 𝑑𝑑 ) corresponds to tendency of the system to form an oblate or a prolate morphology; similarly, the power spectral density profiles of Δ (t) and S( 𝑑𝑑 ) correspond to the vibrational modes of the network and their preferred shearing orientation.
2. Gelation of actin filaments into clusters governed by actin-binding proteins
A complementary approach to study actomyosin networks, which has been used in previous work [10], is to explore the sol-gel phase transition with respect to the molecular connectivity of ABPs and filaments. The 𝑖𝑖 -th cluster at time 𝑑𝑑 is a gelated group of filaments connected by either motors, linkers, or both; we denote as 𝑁𝑁 𝑐𝑐𝑖𝑖 ( 𝑑𝑑 ) the number of filaments in the 𝑖𝑖 -th cluster. By convention, in our framework, a pair of filaments forms the smallest cluster of size two, hence, 𝑁𝑁 𝑐𝑐𝑖𝑖 ( 𝑑𝑑 ) > 1 always holds. We denote 𝑁𝑁 𝑐𝑐 ( 𝑑𝑑 ) , as the total number of clusters in the system at time 𝑑𝑑 , then, we let 𝑁𝑁 gel ( 𝑑𝑑 ) count the total number of gelated filaments in the system: 𝑁𝑁 gel ( 𝑑𝑑 ) = ∑ 𝑁𝑁 𝑐𝑐 𝑖𝑖 ( 𝑑𝑑 ) 𝑁𝑁 𝑐𝑐 ( 𝑠𝑠 ) 𝑖𝑖=1 . ( ) Then, the gelation ratio order parameter reads:
𝑁𝑁� gel ( 𝑑𝑑 ) = filaments 𝑁𝑁 gel ( 𝑑𝑑 ) . ( ) Defining the largest cluster size as: 𝑁𝑁 max ( 𝑑𝑑 ) = max 𝑐𝑐 ( 𝑠𝑠 ) �𝑁𝑁 𝑐𝑐 𝑖𝑖 ( 𝑑𝑑 ) � , ( ) we let the normalized largest cluster size order parameter be: 𝑁𝑁� max ( 𝑑𝑑 ) = filaments 𝑁𝑁 max ( 𝑑𝑑 ) . ( ) The average cluster size reads: 𝜇𝜇 c ( 𝑑𝑑 ) = c ( 𝑠𝑠 ) ∑ 𝑁𝑁 𝑐𝑐 𝑖𝑖 ( 𝑑𝑑 ) 𝑁𝑁 𝑐𝑐 ( 𝑠𝑠 ) 𝑖𝑖=1 , ( ) and the normalized cluster size order parameter is given by: 𝜇𝜇̅ c ( 𝑑𝑑 ) = filaments 𝜇𝜇 c ( 𝑑𝑑 ) . ( ) The standard deviation cluster size is defined as: 𝜎𝜎 c ( 𝑑𝑑 ) = � c ( 𝑠𝑠 ) ∑ �𝑁𝑁 𝑐𝑐 𝑖𝑖 ( 𝑑𝑑 ) − 𝜇𝜇 𝑐𝑐 ( 𝑑𝑑 ) � 𝑐𝑐 ( 𝑠𝑠 ) 𝑖𝑖=1 , ( ) therefore, the normalized standard deviation cluster size order parameter is: 𝜎𝜎� c ( 𝑑𝑑 ) = filaments 𝜎𝜎 c ( 𝑑𝑑 ) . ( ) The reasoning behind dividing gelation-related quantities by 𝑁𝑁 filaments is to map them on the segment [0,1] .
3. Network theory order parameters computed on the graph of filaments The third approach to quantify the morphology of an actomyosin network is to consider the network generated by positions of its filaments. This framework is highly relevant to compare biophysical theories with experiments in which only F-actin is visible via fluorescence microscopy. For this purpose, we harness ideas and order parameters from graph theory or network theory [57, 58] [33]. As mentioned above, the actin graph at time 𝑑𝑑 is 𝐆𝐆 ( 𝑑𝑑 ) ≝ (V(t), E(t)) , given a cutoff distance 𝑑𝑑 cutoff = 200 nm . For higher readability, we occasionally omit the time argument 𝑑𝑑 . Because of the applicability of graph theory to a wide range of scientific fields, there are an enormous number of topological order parameters available in the literature. Common examples include global order parameters such as community structure [47, 59], and clustering coefficient [48], as well as order parameters defined on the level of individual nodes such as centrality [51, 60]. A simple illustration of how the terminology of graph and network theory drawing quantifiable insights on network topology is elaborated in Fig. 2 Although the mechanical elastic and stress characteristics of actomyosin networks have already been studied to a large extent [2, 61-64], embracing the terminology of graph and network theories opens the possibilities to develop a graph-based renormalization group as done for neural networks [65], and potentially develop a theory describing the nonequilibrium theory on the coarse-grained graph representation of complex biological structures such as actomyosin networks [65-68]. FIG. 2. Fundamental graph properties illustrated on toy networks. Two simplified networks are depicted in (a) and (b). Each of the networks comprises five nodes V = {A, B, C, D, E} . The difference between the network in (a) and the network in (b) is the set of edges connects the five nodes in V . While in graph (a) all nodes are connected, in graph (b) there are two connected components {A, C} and {B, D, E} . In (a), the clique number, which is the largest subset of nodes that are all directly connected, is , and the largest clique is {A, B, C} . In contrast, the largest clique in graph (b) is {A, C} . The degree of a node is the number of edges connected to it. For example, in (a), the degree of node B equals . The shortest path between two nodes is the minimal number of edges one needs to hop over to reach from one node to the other. In graph (a), the shortest path between node D and node A is D ↔ E ↔ B ↔ A . When nodes are disconnected as in the case of nodes D and A in (b), the shortest path between them is infinity, by definition. The betweenness of a node is the number of the shortest paths between any other nodes that go through this node. For example, in (a), the node B has the largest betweenness equals , as the shortest paths between A and D, A and E, C and D, and C and E go through B. Community structure in a network indicates the presence of heterogeneity in the density of links in a network: if some nodes are more likely connected with each other than with the rest of the network, they may be said to form a community. This heuristic definition has spawned an enormous literature on the subject of community structure in networks [59]. In this paper, we focus solely on a commonly used method: the Clauset-Newman-Moore greedy modularity maximization method [47], we find the number of communities in 𝐆𝐆 at a specific time point 𝑑𝑑 . The community of a node 𝑣𝑣 ∈ V is noted 𝐶𝐶 𝑣𝑣 ; first, the greedy modularity maximization algorithm assigns a unique community 𝑠𝑠 𝑣𝑣 to each node 𝑣𝑣 ∈ V ; then, the algorithm joins pairs of communities to maximize the modularity: 𝑄𝑄 = | E | ∑ �𝐀𝐀 𝑣𝑣𝑣𝑣 − 𝑘𝑘 𝑣𝑣 𝑘𝑘 𝑤𝑤 | E | � 𝛿𝛿 𝑠𝑠 𝑣𝑣 , 𝑠𝑠 𝑤𝑤 ( 𝑣𝑣 , 𝑣𝑣 ) ∈ E , ( ) where |E| is the number of edges between node 𝑣𝑣 and node 𝑤𝑤 : ( 𝑣𝑣 , 𝑤𝑤 ) ∈ 𝐸𝐸 . 𝐀𝐀 is a symmetric adjacency matrix of 𝐆𝐆 : 𝐀𝐀 𝑣𝑣𝑣𝑣 = �
1, ( 𝑣𝑣 , 𝑤𝑤 ) ∈ 𝐸𝐸 ⋀ ( 𝑤𝑤 , 𝑣𝑣 ) ∈ 𝐸𝐸
0, otherwise , ( ) 𝑘𝑘 𝑣𝑣 is the degree of node 𝑣𝑣 and equals the number of edges connected to 𝑣𝑣 . 𝛿𝛿 𝑎𝑎 , 𝑏𝑏 is Kronecker delta and equals if 𝑎𝑎 = 𝑏𝑏 and otherwise. Using these definitions and the algorithm, we define the order parameter 𝑁𝑁 communites ( 𝑑𝑑 ) to be the number of communities in graph 𝐆𝐆 . In many networks, it has been found that triplets of nodes share edges (form a `triangle’, e.g. the triplet {A, B, C} in Fig. 2) more often than would be expected by random chance [48]. The average clustering coefficient of the network 𝐆𝐆 = (V, E) is: Γ = | V | ∑ 𝛾𝛾 𝑣𝑣𝑣𝑣∈V , ( ) where 𝛾𝛾 𝑣𝑣 is the clustering coefficient of node 𝑣𝑣 : 𝛾𝛾 𝑣𝑣 = ∑ 𝐀𝐀 𝑣𝑣𝑖𝑖 𝐀𝐀 𝑖𝑖𝑖𝑖 𝐀𝐀 𝑣𝑣𝑖𝑖𝑖𝑖 , 𝑖𝑖 𝑘𝑘 𝑣𝑣 ( 𝑘𝑘 𝑣𝑣 −1 ) . ( ) Therefore, Γ ( 𝑑𝑑 ) is the average clustering coefficient order parameter at time 𝑑𝑑 . Note that the clustering coefficient in network theory is related to the density of triangles of edges, and is conceptually distinct from the clustering discussed in the gelation literature. Node centrality has been useful in better understanding the structure of complex networks on the level of individual nodes [51, 60]. There are many possible measures of centrality found in the literature, each measuring the ‘importance’ of a node in a slightly different way. We let 𝑣𝑣 ∈ V be a node in the Θ 𝑣𝑣 ⊆ V connected component of 𝐆𝐆 . Assuming Θ 𝑣𝑣 has 𝑛𝑛 nodes connected with a path between all the nodes, the Wasserman and Faust (WF) closeness centrality [50] of the node 𝑣𝑣 ∈ Θ 𝑣𝑣 ⊆ V is given by: 𝐶𝐶 𝑊𝑊𝑊𝑊 ( 𝑣𝑣 ) = ( 𝑛𝑛−1 ) | V | −1 �∑ 𝑑𝑑 ( 𝑣𝑣 , 𝑢𝑢 ) 𝑢𝑢∈Θ 𝑣𝑣 � −1 , ( ) where 𝑛𝑛 = | Θ 𝑣𝑣 | is the number of nodes in the connected component and 𝑑𝑑 ( 𝑖𝑖 , 𝑗𝑗 ) is the shortest path between the two nodes: 𝑖𝑖 , 𝑗𝑗 ∈ Θ 𝑣𝑣 . The WF centrality is higher for nodes that have short paths to other nodes in the network. Then the mean WF closeness centrality order parameter at time 𝑑𝑑 reads: Ω F ( 𝑑𝑑 ) = ∑ 𝐶𝐶 𝑊𝑊𝑊𝑊 ( 𝑣𝑣 ) 𝑣𝑣∈V . ( ) The eigenvector centrality of node 𝑖𝑖 is the 𝑖𝑖 -th element of the eigenvector 𝐯𝐯 defined by: 𝐀𝐀𝐯𝐯 = 𝜆𝜆 max 𝐯𝐯 , ( ) where 𝜆𝜆 max ≥ 𝜆𝜆 ≥ ⋯ ≥ 𝜆𝜆 | V | is the maximal eigenvalue of the symmetric adjacency matrix 𝐀𝐀 of the graph 𝐆𝐆 . Eigenvector centrality tends to be higher for nodes that are connected to other nodes with high centrality. Then, the average eigenvector centrality is simply: E 𝜆𝜆 max ( 𝑑𝑑 ) = | V ( t )| ∑ 𝜆𝜆 𝑖𝑖1≤𝑖𝑖≤ | V ( t )| . ( ) The average neighbor degree of a node 𝑣𝑣 is evaluated on 𝑣𝑣 ’s neighborhood, the set of the nearest neighbors of 𝑣𝑣 , 𝒩𝒩 ( 𝑣𝑣 ) = { 𝑢𝑢 | ( 𝑢𝑢 , 𝑣𝑣 ) ∈ E} : 𝓀𝓀 𝑛𝑛𝑛𝑛 , 𝑖𝑖 = | 𝒩𝒩 ( 𝑣𝑣 )| ∑ 𝑘𝑘 𝑢𝑢𝑢𝑢∈𝒩𝒩 ( 𝑣𝑣 ) , ( ) where 𝑘𝑘 𝑢𝑢 is the degree of the node 𝑢𝑢 . Like eigenvector centrality, neighbor degree centrality is higher for nodes that are connected to other central nodes, but unlike eigenvector centrality longer paths are not explicitly accounted for. Then we define an order parameter as the mean average neighbor degree: 𝑘𝑘 𝑛𝑛𝑛𝑛 ( 𝑑𝑑 ) = | V | ∑ 𝓀𝓀 𝑛𝑛𝑛𝑛 , 𝑣𝑣𝑣𝑣∈V . ( ) The betweenness centrality of a node (referred to as betweenness) measures how central a node is by how often it acts as a `bridge’ between other pairs of nodes (sketched in Fig. 2). This is quantified by counting the number of shortest paths in which the node participates: 𝛽𝛽 ( 𝑣𝑣 ∈ V) = ∑ 𝑔𝑔 ( 𝑠𝑠 , 𝑠𝑠′ | 𝑣𝑣 ) 𝑔𝑔 ( 𝑠𝑠 , 𝑠𝑠′ ) 𝑠𝑠 , 𝑠𝑠′ ∈V ∗ , ( ) where V ∗ = V ∖ { 𝑣𝑣 } is the set of nodes excluding node 𝑣𝑣 , the function 𝑔𝑔 ( 𝑠𝑠 , 𝑠𝑠′ ) counts the number of geodesics (shortest paths) between node 𝑠𝑠 and node 𝑠𝑠′ , and 𝑔𝑔 ( 𝑠𝑠 , 𝑠𝑠′ | 𝑣𝑣 ) is the number of shortest path between node 𝑠𝑠 and 𝑠𝑠′ such that 𝑣𝑣 lies along those shortest paths. We can also define the degree assortativity coefficient [53], using the definitions from Newman et al. [69, 70] for the nodes’ joint probability matrix. First, we define the degree joint counting matrix 𝐄𝐄 ; this is a symmetric matrix where its 𝐄𝐄 𝑠𝑠𝑙𝑙 element equals the number of edges in the graph that link a vertex of degree 𝑙𝑙 with another vertex of degree 𝑚𝑚 : 𝐄𝐄 𝑠𝑠𝑙𝑙 = |{( 𝑢𝑢 , 𝑣𝑣 ) ∶ 𝑣𝑣 , 𝑢𝑢 ∈ 𝐸𝐸 ⋀ 𝑘𝑘 𝑣𝑣 = 𝑙𝑙 ⋀ 𝑘𝑘 𝑢𝑢 = 𝑚𝑚 }| . ( ) The joint probability matrix is achieved by normalizing 𝐄𝐄 𝑠𝑠𝑙𝑙 with the sum of its elements ∑ 𝐄𝐄 𝑠𝑠𝑙𝑙𝑠𝑠 , 𝑙𝑙 . The degree assortativity coefficient order parameter is given by: 𝜌𝜌 ( 𝑑𝑑 ) = ∑ 𝑠𝑠⋅𝑙𝑙⋅ ( 𝐄𝐄 𝑠𝑠𝑙𝑙 ( 𝑠𝑠 ) −𝑝𝑝 𝑠𝑠 ( 𝑠𝑠 ) 𝑝𝑝 𝑙𝑙 ( 𝑠𝑠 )) 𝑠𝑠 , 𝑙𝑙 𝜎𝜎 𝑝𝑝2 ( 𝑠𝑠 ) , ( ) Where 𝑝𝑝 𝑖𝑖 ≝ ∑ 𝐄𝐄 𝑖𝑖𝑗𝑗𝑗𝑗 = ∑ 𝐄𝐄 𝑗𝑗𝑖𝑖𝑗𝑗 for any possible node’s degree 𝑖𝑖 , and 𝜎𝜎 𝑝𝑝 is the standard variation of the set of the 𝑝𝑝 𝑖𝑖 values, { 𝑝𝑝 𝑖𝑖 } 𝑖𝑖=1 | V | . In addition to the detailed order parameters above, the global details of the network can be considered using coarse measures as well. We consider the graph density, defined by: 𝑑𝑑 𝐆𝐆 ( 𝑑𝑑 ) = | E || V | ⋅ | V−1 | , ( ) which measures the number of edges observed in the network compared to the maximum number that could be observed. We also consider the maximal clique size (sketched in Fig. 2), which is the largest subset of vertices from V such that each two nodes in the subset have an edge between them. F. Potential of Mean Force of Actomyosin Networks
The potential of mean force [71] is the natural logarithm of a probability density function (PDF) defined on a degree of freedom in the system. We let 𝑃𝑃 ( 𝜉𝜉 ) be such a probability density function (PDF) of a parameter 𝜉𝜉 ∈ [0, 1] . Since 𝜉𝜉 is an order parameter, by definition, there exists a mapping from the configuration space of the system onto the segment [0,1] . Subsequently, we can write 𝜉𝜉 to be this function 𝜉𝜉�𝐱𝐱 (t) � which can be evaluated at each time point based on the configuration of the system. Along a trajectory, the probability density function of 𝜉𝜉 is given by: 𝑃𝑃 ( 𝜉𝜉 ) = �𝛿𝛿�𝜉𝜉 − 𝜉𝜉 ( 𝐱𝐱 ) �� trajectory ≝ 𝑄𝑄 ( 𝜉𝜉 ) ∑ 𝑄𝑄 ( 𝜉𝜉 ) 𝜉𝜉 , ( ) where 𝑄𝑄 ( 𝜉𝜉 ) is the partition function to find the system in a state at which the order parameter of interest has the value 𝜉𝜉 . The entropy 𝒮𝒮 ({ 𝜉𝜉 }) [72] [73] on the probability density function of 𝜉𝜉 is given by: 𝒮𝒮 ({ 𝜉𝜉 }) = −𝐾𝐾 𝐵𝐵 ∑ 𝑃𝑃 ( 𝜉𝜉 ′ ) 𝜉𝜉 ′ ln 𝑃𝑃 ( 𝜉𝜉 ′ ) = 𝑘𝑘 𝐵𝐵 ∑ 𝑃𝑃 ( 𝜉𝜉 ′ ) 𝜉𝜉 ′ I( 𝜉𝜉 ′ ) , ( ) where I( 𝜉𝜉 ) ≝ −𝑘𝑘 𝐵𝐵 ln 𝑃𝑃 ( 𝜉𝜉 ) is the Shannon information (self-information) as a function of 𝜉𝜉 . If we let 𝑘𝑘 𝐵𝐵 be unity, both the entropy and the information measure have the nat (natural unit of information) units [74]. The Shannon information function is defined: I( 𝜉𝜉 ) ≝ − 𝑘𝑘 𝐵𝐵 ln 𝑃𝑃 ( 𝜉𝜉 ) ≝ ln � ∑ 𝑄𝑄�𝜉𝜉 ′ � 𝜉𝜉′ 𝑄𝑄 ( 𝜉𝜉 ) � . ( ) We used the python KernelDensity estimator the Scikit Python library [75] to evaluate the self-information I( 𝜉𝜉 ) via the negative logarithm of the PDF 𝑃𝑃 ( 𝜉𝜉 ) . The self-information is also related to the potential of mean force (PMF), noted 𝒰𝒰 ( 𝑇𝑇 , 𝜉𝜉 ) [71, 76] for which the following holds by definition: 𝒰𝒰 ( 𝑇𝑇 , 𝜉𝜉 ) ≝ 𝑇𝑇 I( 𝜉𝜉 ) ≝ − 𝑘𝑘 𝐵𝐵 𝑇𝑇 𝑙𝑙𝑛𝑛 𝑃𝑃 ( 𝜉𝜉 ) ≝ 𝑘𝑘 𝐵𝐵 𝑇𝑇 ln � ∑ 𝑄𝑄�𝜉𝜉 ′ � 𝜉𝜉′ 𝑄𝑄 ( 𝜉𝜉 ) � . ( ) III. RESULTS A. Multilinkers with a high valency enrich arborization of actin bundles
We have conducted the mesoscopic simulations using Cytosim in total of 600 seconds to investigate the role of multivalent linkers on the structural changes in actomyosin. In the simulations, the size of the simulation box is μ m with 1:1:1:1 concentration ratio among filaments, motors, multilinkers, and α -actinin crosslinkers. We have chosen three types of multilinkers with valency of 2, 3, and 6 and provided key snapshots in Figure 3. In the systems with bivalent multilinkers we observed the formation of bundles [77, 78] with a smaller variation in the cluster sizes compared with systems with multilinkers of valency 𝜈𝜈 ≥ as time evolves. In the systems with trivalent multilinkers, we observed nascent arborization of bundles (the appearance of branched actin bundles) at 𝑑𝑑 = 50 s in Fig. 3(a). This phenomenon becomes more prominent in systems with multilinkers of higher valencies as shown in Fig. 3(b) at time 𝑑𝑑 = 600 sec , for a system with hexavalent ( 𝜈𝜈 = 6 ) multilinkers. FIG. 3. Snapshots of simulations from Cytosim of three μ m box systems with 1:1:1:1 concentration ratio of filaments, motors, multilinkers, and α -actinin like crosslinkers. The three systems differ only in the valency of a multilinker: On the most left column in (a) and (b), the multilinkers are bivalent ( 𝜈𝜈 = 2 ). On the middle column the multilinkers are trivalent ( 𝜈𝜈 = 3 ). On the rightest column the multilinkers are hexavalent ( 𝜈𝜈 = 6 ). The time evolves from the top panels to the bottom panels. In (a) the snapshots are in increment of 50 seconds, starting at time 𝑑𝑑 = 𝑑𝑑 , where 𝑑𝑑 = 1 sec . In (b) the snapshots are in increment of 200 seconds, ending at 𝑑𝑑 = 600 sec . B. Correlation of the order parameters from polymer physics, gelation theory, and network theory The rich actomyosin structures shown in Fig 3 have prompted us to seek appropriate order parameters to characterize the complexity and heterogeneity of actomyosin dynamics in space over time. We have compared the order parameters inspired by three domains: polymer physics, gelation theory, and network theory, as described in the Method section. Fig. 4 shows the Pearson correlation of all order parameters evaluated for each actomyosin configuration from the simulations using Cytosim. Each group of order parameters presents a different insight into the actomyosin morphology. Interestingly, the three polymer physics order parameters: radius of gyration ( 𝑅𝑅 𝑔𝑔 ), the shape parameter (S) , and the asphericity ( Δ ) are highly correlated, as all three are related to eigenvalues of the inertia tensor in Eq. (3). 𝑅𝑅 𝑔𝑔 correlates less with S and Δ than S correlates with Δ . These three parameters are nearly uncorrelated to the order parameters from gelation theory or from graph theory. Interestingly, the order parameters from network theory reveal rich information about the connectedness in an actomyosin network. The number of gelated clusters and the node degree assortativity have negative correlation with most order parameters from the gelation and graph theory. The order parameters based on network theory are also clustered together, while the node degree assortativity has a negative correlation with the other gelation and network based order parameters. FIG. 4. Heatmap of the Pearson correlation between the order parameters used for this study. The order parameters are grouped according to the pairwise correlation distance and annotated by their physical order parameter group: polymer physics, gelation theory, or network theory. C. Motor activity regulates temporal changes in the global structure and the shape of actomyosin networks
The fluctuation in the distribution of filament vertices in space is a proxy for the temporal dynamics of the actomyosin networks since the filaments are the scaffold of the network. Using the power spectral density functions (PSD) [56] of the three polymer physics order parameters 𝑅𝑅 𝑔𝑔 , Δ , and S defined in Eq. (3-5): are compared for high and low motor-to-filament ratios in Fig. 5. Each of these three order parameters measures a distinctive mesoscopic reshaping mode of the actomyosin network. PSD of the shape parameter, S, is a mesoscopic assessment of filament sliding modes (Fig 5 (a-b)). The PSD of the asphericity ( Δ ) is a proxy for the rotational modes of filaments (Fig 5(c-d)). The PSD of the radius of gyration ( 𝑅𝑅 𝑔𝑔 ) order parameter corresponds to mesoscopic expansion-contraction modes (Fig 5(e-f)), The filament sliding mode in Fig. 5 (a-b) is more dominant than the rotational mode in (Fig. 5(c-d) because The PSD of the former is one order of magnitude ( ) higher than the PSD of the latter. When comparing the PSD at the two low and high conditions of motor-filament-crosslinker-multilinker concentrations, (1:100:100:100) vs. (1:1:1:1), all the three PSD in the case of high motor-filament ratio is two order of magnitudes ( ) greater than those of low motor-filament ratio. FIG. 5. Power spectral analysis of mesoscopic polymer order parameters. Low and high motor-filament-crosslinker-multilinker ratios, (1:100:100:100) and (1:1:1:1) are examined in systems with multilinkers of various valencies 𝜈𝜈 ∈ {2,3,4,5,6,7} . (a), (c) and (e) show the PSD from the systems with low motor-filament ratios, while (b), (d), and (f) show the PSD from the systems with high motor-filament ratios. We compared the PSD for the shape parameter (S), asphercity ( Δ ) , and the radius of gyration, 𝑅𝑅 𝑔𝑔 . D. Gelation of actomyosin network by actin-binding proteins
Although the shape and structure parameters ( 𝑅𝑅 𝑔𝑔 , S , and Δ ) from polymer physics probe the fluctuations from the configurations of actin filaments in a network, they cannot, however, directly capture the individual characteristics of a gelated, or percolated cluster, in the actomyosin networks [79]. The order parameters from gelation theory, which consider the connectivity between actin filaments through actin-binding proteins, can be used to monitor the development of a percolated cluster. Fig. 6(a-b) shows the normalized largest cluster size over time, 𝑁𝑁� max ( 𝑑𝑑 ), from Eq. (10). It demonstrates that multilinkers with higher valency promote the larger clusters regardless of motor concentration. Fig. 6(c-d) shows the normalized standard deviation of cluster size, 𝜎𝜎� c ( 𝑑𝑑 ) defined in Eq. (14). 𝜎𝜎� c ( 𝑑𝑑 ) increases in the presence of multilinkers with higher valencies at both high and low motor to filament concentration ratios. Fig 6 demonstrates how motors regulate the total number of clusters independent of the valency of multilinkers. The higher the valency of the multilinkers, the faster the multilinkers reduce the number of clusters, 𝑁𝑁 𝑐𝑐 ( 𝑑𝑑 ) in the system at a low level of motors. While at a high level of motors in Fig 6(f), motor activity alone dominates the number of clusters. Fig. 6(g) shows the cumulative distribution functions (CDF) of the average cluster size, 𝜇𝜇̅ c at high (1:1:1:1) and low (1:100:100:100) motor-filament-crosslinker-multilinker concentration ratios in the presence of multilinkers with varying valencies. Multilinkers with valency 𝜈𝜈 = 2 and 𝜈𝜈 = 3 create smaller clusters on average at low motor concentration. However, for a system with multilinkers of 𝜈𝜈 = 4 , both systems with high and low motor-filament ratio have a similar distribution of 𝜇𝜇̅ c . Finally, at 𝜈𝜈 = 5, 𝜈𝜈 = 6, and 𝜈𝜈 = 7 , higher motor activity promotes smaller average cluster sizes. FIG. 6. Gelation of filament clusters formed by crosslinkers, multilinkers, or motors. The condition of "low motor concentration" corresponds to motor-filament-multilinker-crosslinker ratios of 1:100:100:100 as in (a, c, e). The condition of "high motor concentration" corresponds to motor-filament-multilinker-crosslinker ratios of 1:1:1:1 as in (b, d, f). The time evolution of the order parameters (excluding double counting of clusters of the same exact size): the normalized largest cluster as defined in Eq. (10), the time-dependent signal of the order parameter 𝜇𝜇̅ c ( 𝑑𝑑 ) defined in Eq. (12), and the normalized cluster size standard deviation 𝜎𝜎� c ( 𝑑𝑑 ) defined in Eq. (14). (g) The cumulative distribution functions (CDF) of the normalized average cluster size order parameter 𝜇𝜇 c ( 𝑑𝑑 ) , defined in Eq. (11) are plotted in both motor concentration conditions and at the six different multilinker’s valency conditions. Each provide in (a)-(f) is averaged over 30 independent trajectories with the same species parameters with a confidence interval of CI = 90% . E. Constructing graphs/networks using the spatial information of actin filaments
In real world systems, the local microscopic information of actin-binding connectivity to filaments is not always known nor easy to track. In addition, treating actomyosin as percolating network or a gel lacks the information about the hierarchy of connectivity in the scaffold of the network. To further coarse-grain the morphology of actomyosin networks, we introduced the non-local order parameters from network theory where a graph of a network is built by accounting only for the distances between filaments and by overlooking the types of ABPs responsible for connectivity. Hence, we provided an alternative representation of the actomyosin networks in Fig 3 for the two multilinkers of 𝑣𝑣 = 2 and 𝑣𝑣 = 6 at 𝑑𝑑 = 1 s and at 𝑑𝑑 = 600 s in Fig. 7, where the node’s size and its purple color intensity are scaled with its betweenness 𝛽𝛽 ( 𝑣𝑣 ) defined in Eq. (25). The green-white colors correspond to whether two nodes have a connective edge or not. In this graph representation of the actomyosin networks, the nodes are filaments and an edge between two nodes exists if they are within a cutoff distance 𝑑𝑑 cutoff = 200 nm . Fig. 7(a) shows the network of at 𝑑𝑑 =1 sec where the multilinkers are bivalent ( 𝜈𝜈 = 2 ). Fig. 7(b) shows the network layout at 𝑑𝑑 = 1 sec for a system with hexavalent multilinkers ( 𝜈𝜈 = 6 ). Fig. 7(c) and Fig. 7(d) depict the two networks at 𝑑𝑑 = 6 sec . The topologies of the two networks appear distinctive after 600 seconds of evolution. These distinctive features have purposefully captured by the order parameters from graph theory: For the graph representing the bivalent multilinkers ( 𝜈𝜈 = 2 ), the graph density is , the average node degree is , and the average path length 𝑖𝑖𝑠𝑠 whereas for the graph representing the hexavalent multilinkers ( 𝜈𝜈 = 6 ), the graph density of , the average node degree of , and the average path length of . In the following sections, we will explore the distinctive phases of actomyosin networks, that are not apparent in previous sections, using the order parameters from network theory. FIG. 7. Graph representations of actomyosin networks using Gephi 0.92 [80]. 𝐆𝐆 = (V, E) , where the nodes V represent the filaments whose adjacency matrix is binary and an edge between filament 𝑖𝑖 and filament 𝑗𝑗 exists in E if the Euclidean distance between their center of masses is less than 𝑑𝑑 cutoff = 200 nm . The multilinker valency in the left column (a, c) is 𝜈𝜈 = 2 , and that in the right colume (b, d) is 𝜈𝜈 = 6 . The top panels (a, b) are snapshots at t= , and the bottom panels (c, d) are snapshots at time 𝑑𝑑 = 600 s . The ratio of filaments, motors, multilinkers, crosslinkers is 1:1:1:1 in all panels. We produced these graphs by following a visualization procedure: First we produced a random network layout with OpenOrd [81] and utilized Yifan Hu’s algorithm [82] with the default parameters in Gephi 0.92. The betweenness of a node counts the number of times it lies on the shortest path between any two other nodes in Eq. (24). F. Analysis of actin filament networks viewed by order parameters from graph theory
The topology of actomyosin network becomes increasingly complex with time, as shown in Fig. 7(c-d). We are interested in characterizing its ramifications of neighboring filaments. Network theory provides a layperson’s description on the ranks of the social hierarchy in a network built by these connected filaments. In this sense, clique numbers, number of communities, and clustering numbers report the status of a node, or filament, in a growing dendritic network. Graph theory allows us to view the network from a fine-grained (zoomed-in) to a coarse-grained (zoomed-out) perspective. The strength of network theory to describe scales beyond the scope of ABPs provide the social status of a node (or filament) in a hierarchical actomyosin network. The interplay of multilinker valency and motor proteins on the morphology of an actomyosin network with several order parameters such as clique number, degree assortativity, and graph density in Fig 8, 9 and 10, respectively. Lastly, we exploited one of these order parameters in the graph theory to characterize the complex energy landscapes in terms of potential of mean forces in statistical physics in Fig 11. We first zoomed in an actomyosin network at the fine-grained level where multivalency differentiates topologies at either high or low motor content in Fig 8. At low motor content, filaments produce a similar kind of bundles at all valencies from 𝜈𝜈 = 2 to 𝜈𝜈 = 7 according to Fig. 8(a, c, e). Whereas, the differences in network topology among multilinkers with various valencies is most prominent at high motor activity as the average clustering profile (Fig. 8b), the clique number (Fig. 8d), and the number of communities (Fig. 8f) diverse over time. Prominently, the number of communities [47] that maximize the modularity in Eq. (15) is the highest for 𝜈𝜈 = 7 in Fig 8f. For systems with multilinker of higher valencies and this is due to the effect of high valency on the average degree of the nodes in the graph as shown in Fig. 7(c-d) which separate the cluster assignment of filaments into separated clusters so the modularity in Eq. (15) remains positive and close to zero. When we zoomed out the network at a coarse-grained level to examine the node degree coefficient that depicts hierarchical status of a node (filament) in a network. Higher valency multilinkers reduce the node degree coefficient defined in Eq. (27) both at 1:100:100:100 and 1:1:1:1 motor-filament-crosslinker-multilinker ratios as shown in Fig. 9(a) and Fig. 9(b), respectively. Because multilinkers with multivalency 𝜈𝜈 ≥ , (Fig. 9(b)) at valencies 𝜈𝜈 ≥ , weighed more important socially in a network, motors in high concentration (Fig. 9(c-d)) through these highly ranked multilinkers effectively reduces the assortativity and causes arborization of actin bundles. Indeed, motors reduce the assortativity of the network from to for the cases of hexavalent ( 𝜈𝜈 = 6 ) and heptavalent ( 𝜈𝜈 = 7 ) multilinkers. G. Building Network/graph theory analysis of actin filament networks
Note that the network/graph between filaments is built by accounting only for the distances between filaments and neglecting the ABPs responsible for connectivity. Fig. 7 shows four networks where the nodes are filaments and an edge between two nodes exists if they are within a cutoff distance 𝑑𝑑 cutoff = 200 nm . Fig. 7(a) shows the network of at 𝑑𝑑 = 1 sec where the multilinkers are bivalent ( 𝑣𝑣 = 2 ). Fig. 7(b) shows the network layout at 𝑑𝑑 = 1 sec for a system with hexavalent multilinkers ( 𝜈𝜈 = 6 ). Fig. 7(c) and Fig. 7(d) depict the networks at 𝑑𝑑 = 600 sec . The topologies of the two networks look different after 600 seconds of evolution; the network with bivalent multilinkers has a graph density of , an average node degree of , and an average path length of compared with graph density of , an average node degree of , and an average path length of in the system with hexavalent multilinkers. At 1:100 motor-filament concentration ratio together with multilinkers with valency 𝜈𝜈 ≥ , the network of the filaments produces bundles similar to systems with bivalent multilinkers, as seen in Fig. 8(a) and Fig. 8(e). The differences in network topology between the different types of multilinkers is most prominent at high motor activity as shown for the average clustering profile in Fig. 8(b), the clique number in Fig. 8(d), and number of communities in Fig. 8(f). The number of communities [47] that maximizes the modularity in Eq. (15) is higher for systems with multilinker of higher valencies, This is due to the effect of high valency on the average degree of the nodes in the graph as shown in Fig. 7(c-d) which separate the cluster assignment of filaments into distinct clusters. Consequently, the modularity in Eq. (15) remains positive and close to zero. Higher valency multilinkers reduce the node degree coefficient defined in Eq. (27) both at 1:100 and 1:1 motor-filament ratios as shown in Fig. 9(a) and Fig. 9(b), respectively. As shown in Fig. 9(b), at valencies 𝜈𝜈 ≥ , high motor concentration enhances the assortativity reduction and causes arborization of actin bundles; motors reduce the assortativity of the network from to for the cases of hexavalent ( 𝜈𝜈 = 6 ) and heptavalent ( 𝜈𝜈 = 7 ) multilinkers. FIG. 8. Order parameters from network theory measuring connectivity in the presence of multilinkers for two systems with 1000 nodes (filaments). The left column (a, c, e) represents the system with (1:100:100:100) concentration ratios among motors, filaments, multilinkers, and cross linkers. The right column (b, d, f) represents the system with (1:1:1:1) concentration ratios among motors, filaments, multilinkers, and cross linkers. The average clustering, defined in Eq. (17) is plotted as function of time in (a, b). The clique number of the network as function of time is plotted in (c, d). The number of communities [47] of the network is plotted as function of time in (e, f). Each profile from (a) to (f) is averaged over 30 independent trajectories with a confidence interval of
CI = 90% . FIG. 9. Node degree assortativity coefficient of an actomyosin network in the presence multilinkers with varying valencies. The node degree assortativity defined in Eq. (27) is plotted as a function of time at 1:100:100:100 concentration ratios between motors, filaments, multilinkers, and crosslinkers in (a) and at 1:1:1:1 ratios in (c). Each profile in (a) and (c) is averaged pver 30 independent trajectories with a confidence interval of
CI = 90% . The snapshots in (b) and (d) are of μ m cubic systems with multilinkers of valency 𝜈𝜈 = 6 at low and high motor content at time 𝑑𝑑 = 600 𝑠𝑠 , respectively. H. Network theory order parameters reveal rich topologies in a dendritic actomyosin network
FIG. 10. Graph density of actomyosin networks in the presence of multilinkers with varying valencies. The graph density defined in Eq. (28) is plotted as a function of time at 1:100:100:100 concentration ratios between motors, filaments, multilinkers, and crosslinkers in (a) and at 1:1:1:1 ratios in (b). Each profile in (a) and (b) is averaged of 30 independent trajectories with the same species parameters with a confidence interval of
CI = 90% . Beyond the order parameters counting the node (or filament) connectivity, the non-local order parameter in graph theory such as graph density profile reveals the spatial heterogeneity of actomyosin networks at a coarse-grained level, which tracks the branching dendrites in an overall network. In Fig 10(a), graph density shows that even at low motor activity, the presence of multilinkers with 𝜈𝜈 > 2 increases the density by 3-fold. At high motor content in Fig. 10(b), suggests the notion that the dynamics of the network has more than two phases as in the case of gelation and solution. For the case of multilinker with 𝜈𝜈 = 3 , the graph density shows that there are two stages of evolution characterized by a double sigmoidal curve: the graph density is less than 0.1 before 𝑑𝑑 = 200 s and abruptly grows to 0.2 after t>300. This double sigmoidal curve exposes that there more than two phases, gelation and solution, in the dendritic growth of actomyosin networks. Therefore, tracking only percolation using the order parameters in gelation theory may be inadequate. We report that using a coarse-grained level of order parameter in graph theory is highly sufficient in revealing the hierarchical phases in a dendritic growth of actomyosin networks. Because these order parameters capture the hierarchy of actin filaments in an actomyosin network, they are useful in revealing the physics of dendritic growth over time. It is desirable to develop them into the order parameters that describe the complexity of an energy landscape. To show such a possibility, we show the potential of the mean force (PMF) defined in Eq. (32) as a function of several of these order parameters in Fig 11. for the sizes of gel clusters and compares it with the PMF of the average node neighbor degree of each node among a trajectory of
600 s . Fig. 11 does support the observation that with network theory (Fig. 11(a)) order parameters there are more than two phases (noted by the number of local minima) of actomyosin with multilinkers compared with a two state of sol-gel profile in Fig. 11(b). FIG. 11. Potential of Mean Force (PMF) of (a) network and (b) gelation theory parameters. The potential of mean force) as defined in Eq. (32) of the probability density function over the cluster sizes in (a) of systems with multilinkers of different valencies. The PMF of the probability density function of the average neighbor degree of each node is plotted in (b). Both (a) and (b) were computed at 1:100:100:100 concentration ratios between motors, filaments, multilinkers, and crosslinkers. I. State diagrams of actomyosin networks in the plane of order parameters from network theory
With the proper order parameters connecting local and non-local features in a complex network established from network theory, we have outlined the state diagrams in Fig 12 that allow further experimental validations on our hypothesis that both motor activity and the valency of the multilinker in the system regulate the dendritic state of the actomyosin network. The state diagrams in Fig. 12(a-b) are computed by averaging over all the last 20 seconds of simulations for which a specific ratio of motor-filament and multilinker-filament holds. A point in Fig. 12(c) is computed by averaging over the last 20 seconds of simulations for which the multilinker valency is fixed, and multilinker-filament ratio is fixed, when a gaussian filter is applied along the discrete valency axis. First, we show that the state diagram of the degree of assortativity at 2:1 motor-filament concentration ratio or above the degree assortativity is reduced as shown in Fig. 12(a). At high multilinker-filament concentration ratio above 2:1 the degree assortativity reduces as well. Both high concentration of multilinkers and high activity of motors are necessary to regulate the correlation between nodes in the network. At multilinker-filament concentration ratio above 4:1 as shown in Fig. 12(a) the motor activity does not attenuate the assortativity as in the case where the motor-filament concentration ratio is between 1:2 and 2:1. Fig. 12(b) shows that the graph density of the network does not necessarily increase by adding more multilinkers, nor by adding more motors. Fig. 12(c) shows that at 1:1 multilinker-filament concentration ratios and lower, the clique number increases with multilinker with valency 𝜈𝜈 ≥ as also shown for the 1:1 motor-filament concentration ratio in Fig. 8(d). FIG. 12. State diagrams of order parameters from network theory. The average over the steady state values of the node degree assortativity coefficient, defined in Eq. (27) are in (a). The average over the steady state values of the graph density, defined in Eq. (28), are plotted in (b). The average of the steady state values over the clique number normalized by the number of nodes (filament) in the network is shown in (c). In order to create a continuous plot in (c) along the Valency of Multilinker axis a gaussian filter was used to smooth the signal along the valency of multilinker (x-axis). IV. DISCUSSION AND CONCLUSION A. Network theory approach is necessary to characterize the dendritic features on actomyosin networks
We have simulated the growth of actomyosin networks mediated by actin binding proteins using Cytosim. We have utilized several order parameters from three domains of sciences to characterize the dendritic growth of actomyosin. We found that order parameters from polymer physics as well as gelation theory are valuable approaches, but they lack the ability to characterize the topology of a cluster. Using network theory allowed us to see learn more about the complicated dendritic dynamics of actomyosin. The characterization of a topology in a network requires fine-grained descriptions that reveal the connectivity between nodes (or filaments) as well as coarse-grained descriptions that show the overall morphology (scaffold) of hierarchical bundles. In a fine-grained level, the clique number for example from network theory shows the organization of filaments by the valency of the multilinkers and by the motor activity. In a coarse-grained level, the graph density from network theory allowed us to predict the effects of the valency multilinkers versus motor concentration on the topology of the networks through a state diagram. Such descriptions are naturally across multiple length and time scales, and helpful in connecting microscopic properties to macroscopic phenomena that enable experimental validation. For example, when the content of motors or multilinkers increases, at a fine-grained level the order parameters such as the clique number of the network as well as the degree correlation between nodes reduce. At a coarse-grained level, the system becomes less dense by turning hierarchically organized, as the graph density is reduced. This finding resembles the known effect from experimental work that adding more crosslinkers into systems prevents the system from contracting [83]. Our work shows that many interesting features of actomyosin networks are captured by focusing on the nodes (or filaments) from the level of connectivity through actin-binding proteins to a larger scale by knowing only the distances between actin filaments without specific knowledge of actin-binding proteins. Such reduced representation can ease the analysis of the dynamics and structures of the actomyosin networks, reflecting the reorganizations caused to motor and actin-binding proteins. Such a framework to evaluate actomyosin networks is valuable also in the context of analysis of live high-resolution experiments where tracking the location of filaments is usually done. B. Actomyosin structures mediated by multivalent actin-binding proteins shape the synaptic plasticity in a dendritic spine
CaMKII is a multi-valent actin-binding protein that serves a dual role as a structural protein that organizes actin filaments and a signaling protein that mediates calcium signaling. Regarding its structural role, Khan’s as well as Waxham’s group have shown that CaMKII is a multivalent linker that forms actin bundles [84] as well as junctions [5]. CaMKII enriches dendritic actin geometries and augments actin polymerization reshaping the dendritic spine from a filipodia conformation[85] to a mushroom-like structure. CaMKII binds actin in high affinity that maintains the stability of actin bundles [84]over a long period of time in terms of days or months. On the other hand, once calcium-bound calmodulin (CaM) activates CaMKII in a CaMKII-bound actin bundles, CaMKII quickly releases actin in less than a second [29, 86]. A recent work by Wang [4] solved this seemingly conflicting role of CaMKII by showing that calcium-bound CaM significantly weakens the affinity of CaMKII/F-actin complexes by striping one of CaMKII domain away from the actin through binding on the same interfacial residues. The dual role of CaMKII permits it to control the growth as well as to maintain the stability of actin filaments. The unique feature of multivalency in CaMKII has long been speculated to be important in the formation of dendritic filaments and actin bundles. Here, we have run simulations with Cytosim and investigated the impact of multivalency in actomyosin dynamics. With the non-local order parameters such as graph density as well as node degree assortativity from graph theory, we identified the importance of high valency multilinkers such as CaMKII and with high motor content in the dendritic growth of actin filaments which spread in space over time, arborizing from actin bundles to many actin branches. With these useful non-local order parameters from network theory, we have provided state diagrams that permit experimental validations.
C. Concluding Remarks In summary, we have introduced a computational mesoscopic model for multilinkers built into Cytosim. We ran the simulations of actomyosin dynamics in the presence of multilinkers with varying binding valencies to actin filaments. We analyzed the effects of motor activity and the valency of multilinkers on the morphologies of actomyosin with the order parameters from polymer physics, gelation theory and network theory. Our results elucidate the relevance of these order parameters for the study of the structure and dynamics of actomyosin networks. We suggest a possible reduced space from the viewpoint of network theory to characterize complex actomyosin networks. We found that systems that contain multilinkers with valency 𝜈𝜈 ≥ in the presence of high motor concentration exhibit arborized filaments which become apparent in a mathematical network representation. Although this is only a start, we hope that this will help us to understand the role of these molecules in vivo and make experimentally testable predictions. More generally, we expect this work to motivate future directions to embrace network theory as a tool to analyze experiments as well as to establish a quantitative understanding of cytoskeletal systems. We believe that the representation of actomyosin as graphs would be useful beyond the scope of our analysis. The realization that actomyosin could think of as graphs will allow to apply automate the process of coarse graining using graph neural network [87]. Then using nonequilibrium physics theories [68, 88, 89] and statistical physics on probabilistic graphical models [67] could help to develop a theory of nonequilibrium ensembles and landscape reconstruction [40, 41] for actomyosin networks. ACKNOWLEDGMENTS
We thank the members from the Center for Theoretical Biological Physics (CTBP) at Rice University and its memory-focused group. YE and MSC would like to thank Dr. Donald Kouri, Dr. Andre C. Barato, and Dr. Kresimir Josic from the University of Houston for their helpful suggestions and stimulating discussions. Moreover, we thank the Research Computing Data Core (RCDC) Center at the University of Houston for the computational resources. Ultimately, we thank the National Science Foundation for their funding support (CHE:1743392, PHY:1427654, OAC:1531814). Lastly, FJN thanks for his support by the Gatsby Charitable Foundation. The simulations were performed using the Open Source project Cytosim, and the source code is available at https://gitlab.com/f.nedelec/cytosim/-/tags/multilinkers. [1] S. Stam, S.L. Freedman, S. Banerjee, K.L. Weirich, A.R. Dinner, M.L. Gardel, Filament rigidity and connectivity tune the deformation modes of active biopolymer networks, Proceedings of the National Academy of Sciences, 114 (2017) E10037-E10045. [2] M. Malik-Garbi, N. Ierushalmi, S. Jansen, E. Abu-Shah, B.L. Goode, A. Mogilner, K. Keren, Scaling behaviour in steady-state contracting actomyosin networks, Nature Physics, 15 (2019) 509. [3] M.D. Rubio, R. Johnson, C.A. Miller, R.L. Huganir, G. Rumbaugh, Regulation of synapse structure and function by distinct myosin II motors, Journal of Neuroscience, 31 (2011) 1448-1460. [4] Q. Wang, M. Chen, N.P. Schafer, C. Bueno, S.S. Song, A. Hudmon, P.G. Wolynes, M.N. Waxham, M.S. Cheung, Assemblies of calcium/calmodulin-dependent kinase II with actin and their dynamic regulation by calmodulin in dendritic spines, Proceedings of the National Academy of Sciences, DOI (2019) 201911452. [5] S. Khan, K.H. Downing, J.E. Molloy, Architectural Dynamics of CaMKII-Actin Networks, Biophysical journal, 116 (2019) 104-119. [6] T. Nakagawa, J.A. Engler, M. Sheng, The dynamic turnover and functional roles of α -actinin in dendritic spines, Neuropharmacology, 47 (2004) 734-745. [7] T.H. Tan, M. Malik-Garbi, E. Abu-Shah, J. Li, A. Sharma, F.C. MacKintosh, K. Keren, C.F. Schmidt, N. Fakhri, Self-organized stress patterns drive state transitions in actin cortices, Science advances, 4 (2018) eaar2847. [8] J.T. Parsons, A.R. Horwitz, M.A. Schwartz, Cell adhesion: integrating cytoskeletal dynamics and cellular tension, Nat. Rev. Mol. Cell Biol., 11 (2010) 633. [9] J.M. Belmonte, M. Leptin, F. Nédélec, A theory that predicts behaviors of disordered cytoskeletal networks, Molecular systems biology, 13 (2017) 941. [10] J. Liman, C. Bueno, Y. Eliaz, N.P. Schafer, M.N. Waxham, P.G. Wolynes, H. Levine, M.S. Cheung, The role of the Arp2/3 complex in shaping the dynamics and structures of branched actomyosin networks, Proceedings of the National Academy of Sciences, DOI 10.1073/pnas.1922494117(2020) 201922494. [11] B.A. Smith, K. Daugherty-Clarke, B.L. Goode, J. Gelles, Pathway of actin filament branch formation by Arp2/3 complex revealed by single-molecule imaging, Proceedings of the National Academy of Sciences, 110 (2013) 1285-1290. [12] Y.C. Fontela, H. Kadavath, J. Biernat, D. Riedel, E. Mandelkow, M. Zweckstetter, Multivalent cross-linking of actin filaments and microtubules through the microtubule-associated protein Tau, Nature communications, 8 (2017) 1981. [13] B. DuBoff, J. Götz, M.B. Feany, Tau promotes neurodegeneration via DRP1 mislocalization in vivo, Neuron, 75 (2012) 618-632. [14] F. Nedelec, D. Foethke, Collective Langevin dynamics of flexible cytoskeletal fibers, New Journal of Physics, 9 (2007) 427. [15] R. Loughlin, R. Heald, F. Nédélec, A computational model predicts Xenopus meiotic spindle organization, The Journal of cell biology, 191 (2010) 1239-1249. [16] E.A. Nimchinsky, B.L. Sabatini, K. Svoboda, Structure and function of dendritic spines, Annual review of physiology, 64 (2002) 313-353. [17] M. Park, J.M. Salgado, L. Ostroff, T.D. Helton, C.G. Robinson, K.M. Harris, M.D. Ehlers, Plasticity-induced growth of dendritic spines by exocytic trafficking from recycling endosomes, Neuron, 52 (2006) 817-830. [18] Y. Yang, X.-b. Wang, M. Frerking, Q. Zhou, Spine expansion and stabilization associated with long-term potentiation, Journal of Neuroscience, 28 (2008) 5740-5751. [19] K. Zito, G. Knott, G.M. Shepherd, S. Shenolikar, K. Svoboda, Induction of spine growth and synapse formation by regulation of the spine actin cytoskeleton, Neuron, 44 (2004) 321-334. [20] H. Kasai, M. Fukuda, S. Watanabe, A. Hayashi-Takagi, J. Noguchi, Structural dynamics of dendritic spines in memory and cognition, Trends in neurosciences, 33 (2010) 121-129. [21] J. Ryu, L. Liu, T.P. Wong, D.C. Wu, A. Burette, R. Weinberg, Y.T. Wang, M. Sheng, A critical role for myosin IIb in dendritic spine morphology and synaptic function, Neuron, 49 (2006) 175-182. [22] K.U. Bayer, H. Schulman, CaM Kinase: Still Inspiring at 40, Neuron, 103 (2019) 380-394. [23] N.E. Erondu, M.B. Kennedy, Regional distribution of type II Ca2+/calmodulin-dependent protein kinase in rat brain, Journal of Neuroscience, 5 (1985) 3270-3277. [24] S. Khan, I. Conte, T. Carter, K.U. Bayer, J.E. Molloy, Multiple CaMKII binding modes to the actin cytoskeleton revealed by single-molecule imaging, Biophys. J . 111 (2016) 395-408. [25] A.J. Baucum, B.C. Shonesy, K.L. Rose, R.J. Colbran, Quantitative proteomics analysis of CaMKII phosphorylation and the CaMKII interactome in the mouse forebrain, ACS chemical neuroscience, 6 (2015) 615-631. [26] H. Sanabria, M.T. Swulius, S.J. Kolodziej, J. Liu, M.N. Waxham, βCaMKII regulates actin assembly and structure, Journal of Biological Chemistry, 284 (2009) 9770-9780. [27] P.-M. Lledo, G.O. Hjelmstad, S. Mukherji, T.R. Soderling, R.C. Malenka, R.A. Nicoll, Calcium/calmodulin-dependent kinase II and long-term potentiation enhance synaptic transmission by the same mechanism, Proceedings of the National Academy of Sciences, 92 (1995) 11175-11179. [28] K.-I. Okamoto, R. Narayanan, S.H. Lee, K. Murata, Y. Hayashi, The role of CaMKII as an F-actin-bundling protein crucial for maintenance of dendritic spine structure, Proceedings of the National Academy of Sciences, 104 (2007) 6418-6423. [29] J. Lisman, R. Yasuda, S. Raghavachari, Mechanisms of CaMKII action in long-term potentiation, Nature reviews neuroscience, 13 (2012) 169. [30] M. Levitt, A. Warshel, Computer simulation of protein folding, Nature, 253 (1975) 694. [31] M. Vendruscolo, E. Kussell, E. Domany, Recovery of protein structure from contact maps, Folding and Design, 2 (1997) 295-306. [32] S.S. Rao, M.H. Huntley, N.C. Durand, E.K. Stamenova, I.D. Bochkov, J.T. Robinson, A.L. Sanborn, I. Machol, A.D. Omer, E.S. Lander, A 3D map of the human genome at kilobase resolution reveals principles of chromatin looping, Cell, 159 (2014) 1665-1680. [33] R. Albert, A.-L. Barabási, Statistical mechanics of complex networks, Reviews of modern physics, 74 (2002) 47. [34] A.-L. Barabási, Linked: The new science of networks, AAPT, 2003. [35] J.E. Kollmer, K.E. Daniels, Betweenness centrality as predictor for forces in granular packings, Soft matter, 15 (2019) 1793-1798. [36] G. Evenbly, G. Vidal, Tensor network renormalization, Physical review letters, 115 (2015) 180405. [37] Y. Levine, O. Sharir, N. Cohen, A. Shashua, Quantum entanglement in deep learning architectures, Physical review letters, 122 (2019) 065301. [38] M.E. Newman, Spread of epidemic disease on networks, Physical review E, 66 (2002) 016128. [39] M. Kitsak, L.K. Gallos, S. Havlin, F. Liljeros, L. Muchnik, H.E. Stanley, H.A. Makse, Identification of influential spreaders in complex networks, Nature physics, 6 (2010) 888. [40] N. López-Alamilla, M.W. Jack, K. Challis, Reconstructing free-energy landscapes for nonequilibrium periodic potentials, Phys Rev E, 97 (2018) 032419. [41] Y. Fan, T. Iwashita, T. Egami, Energy landscape-driven non-equilibrium evolution of inherent structure in disordered material, Nature communications, 8 (2017) 15417. [42] F. Gittes, B. Mickey, J. Nettleton, J. Howard, Flexural rigidity of microtubules and actin filaments measured from thermal fluctuations in shape, The Journal of cell biology, 120 (1993) 923-934. [43] E.F. Pettersen, T.D. Goddard, C.C. Huang, G.S. Couch, D.M. Greenblatt, E.C. Meng, T.E. Ferrin, UCSF Chimera—a visualization system for exploratory research and analysis, Journal of computational chemistry, 25 (2004) 1605-1612. [44] L.C. Griffith, Regulation of calcium/calmodulin-dependent protein kinase II activation by intramolecular and intermolecular interactions, Journal of Neuroscience, 24 (2004) 8394-8398. [45] R.I. Dima, D. Thirumalai, Asymmetry in the shapes of folded and denatured states of proteins, The Journal of Physical Chemistry B, 108 (2004) 6564-6570. [46] C. Hyeon, R.I. Dima, D. Thirumalai, Size, shape, and flexibility of RNA structures, The Journal of chemical physics, 125 (2006) 194905. [47] A. Clauset, M.E. Newman, C. Moore, Finding community structure in very large networks, Physical review E, 70 (2004) 066111. [48] J. Saramäki, M. Kivelä, J.-P. Onnela, K. Kaski, J. Kertesz, Generalizations of the clustering coefficient to weighted complex networks, Physical Review E, 75 (2007) 027105. [49] S. Wasserman, K. Faust, Social network analysis: Methods and applications, Cambridge university press1994. [50] P. Crucitti, V. Latora, S. Porta, Centrality measures in spatial networks of urban streets, Physical Review E, 73 (2006) 036125. [51] P. Bonacich, Power and centrality: A family of measures, American journal of sociology, 92 (1987) 1170-1182. [52] A. Barrat, M. Barthelemy, R. Pastor-Satorras, A. Vespignani, The architecture of complex weighted networks, Proceedings of the national academy of sciences, 101 (2004) 3747-3752. [53] M.E. Newman, Mixing patterns in networks, Physical Review E, 67 (2003) 026126. [54] A. Hagberg, P. Swart, D. S Chult, Exploring network structure, dynamics, and function using NetworkX, Los Alamos National Lab.(LANL), Los Alamos, NM (United States), 2008. [55] P. Talkner, R.O. Weber, Power spectrum and detrended fluctuation analysis: Application to daily temperatures, Physical Review E, 62 (2000) 150. [56] P. Welch, The use of fast Fourier transform for the estimation of power spectra: a method based on time averaging over short, modified periodograms, IEEE Transactions on audio and electroacoustics, 15 (1967) 70-73. [57] R. Cohen, S. Havlin, Complex networks: structure, robustness and function, Cambridge university press2010. [58] P. Erdős, A. Rényi, On the evol ution of random graphs, Publ. Math. Inst. Hung. Acad. Sci, 5 (1960) 17-60. [59] S. Fortunato, Community detection in graphs, Physics reports, 486 (2010) 75-174. [60] M. Newman, Networks, Oxford university press2018. [61] M. Gardel, J.H. Shin, F. MacKintosh, L. Mahadevan, P. Matsudaira, D. Weitz, Elastic behavior of cross-linked and bundled actin networks, Science, 304 (2004) 1301-1305. [62] F. MacKintosh, J. Käs, P. Janmey, Elasticity of semiflexible biopolymer networks, Physical review letters, 75 (1995) 4425. [63] C. Storm, J.J. Pastore, F.C. MacKintosh, T.C. Lubensky, P.A. Janmey, Nonlinear elasticity in biological gels, Nature, 435 (2005) 191. [64] M.L. Gardel, F. Nakamura, J.H. Hartwig, J.C. Crocker, T.P. Stossel, D.A. Weitz, Prestressed F-actin networks cross-linked by hinged filamins replicate mechanical properties of cells, Proceedings of the National Academy of Sciences, 103 (2006) 1762-1767. [65] S.-H. Li, L. Wang, Neural network renormalization group, Physical review letters, 121 (2018) 260601. [66] G. Gallavotti, E.G.D. Cohen, Dynamical ensembles in nonequilibrium statistical mechanics, Physical review letters, 74 (1995) 2694. [67] D. Koller, N. Friedman, Probabilistic graphical models: principles and techniques, MIT press2009. [68] C. Jarzynski, Nonequilibrium equality for free energy differences, Physical Review Letters, 78 (1997) 2690. [69] D.S. Callaway, J.E. Hopcroft, J.M. Kleinberg, M.E. Newman, S.H. Strogatz, Are randomly grown graphs really random?, Physical Review E, 64 (2001) 041902. [70] M.E. Newman, Assortative mixing in networks, Physical review letters, 89 (2002) 208701. [71] M. Tuckerman, Statistical mechanics: theory and molecular simulation, Oxford university press2010. [72] E.T. Jaynes, Information theory and statistical mechanics, Physical review, 106 (1957) 620. [73] U. Seifert, Stochastic thermodynamics, fluctuation theorems and molecular machines, Reports on progress in physics, 75 (2012) 126001. [74] C.E. Shannon, A mathematical theory of communication, Bell system technical journal, 27 (1948) 379-423. [75] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, Scikit-learn: Machine learning in Python, Journal of machine learning research, 12 (2011) 2825-2830. [76] D. Chandler, Introduction to modern statistical, Mechanics. Oxford University Press, Oxford, UK, DOI (1987). [77] L.G. Tilney, P.S. Connelly, K.A. Vranich, M.K. Shaw, G.M. Guild, Why are two different cross-linkers necessary for actin bundle formation in vivo and what does each cross-link contribute?, The Journal of cell biology, 143 (1998) 121-133. [78] J.D. Winkelman, C. Suarez, G.M. Hocky, A.J. Harker, A.N. Morganthaler, J.R. Christensen, G.A. Voth, J.R. Bartles, D.R. Kovar, Fascin- and αand α