{'fraction_non_alphanumeric': 0.03537081385182651, 'fraction_numerical': 0.01225909453757555, 'mean_word_length': 4.817669172932331, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Multivalent actin-binding proteins augment the variety of morphologies in actomyosin networks

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 Ca2+/Calmodulin-dependent Kinase II (CaMKII). 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 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 2+ /calmodulindependent 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.

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 10 −3 μm 3 and 1 μm 3 [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 2+ /calmodulin-dependent kinase II [22] (CaMKII) holoenzyme protein serves a dual role as a Ca 2+ 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.

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 2 [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: 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 �2 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].

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: where 0 ≥ 0 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 ( 0 ).
The abscissa is increased by = max motor �1 − �, where max motor 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 = 0 exp � | | 0 � to model dissociation from the fiber.
Motors have constant binding on motor and unbinding rates off motor . Namely, they work in the limit of constant unbinding rate ( 0 → ∞). In addition, we set max motor = 0.2 μm s −1 and = 6 pN.
The bivalent crosslinker and multilinkers have the unbinding and binding rates: on linker and off linker .
The interaction of a linker and filaments is assumed to be Hookean with zero resting length and a stiffness of 50 pNμm −1 . , 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 6 nm 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: where is the distance between two interacting elements, 0 is the distance at which the elements are at equilibrium, and 1 is an additional distance that confers an attractive component to the interaction in the range 0 < < 0 + 1 . We use only a repulsive steric effect between multilinkers with 0 = 5 nm, − = 0, and + = 500 pNμm −1 and thus 1 = 0.

F. General Simulation's Setup Parameters
Our system is a

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 (Rg), 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.

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 − 1 4 ≤ S ≤ 2, and the asphericity measure 0 ≤ Δ ≤ 1. 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 : 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, 1 ≥ 2 ≥ 3 of at time : We let ̅ ( ) = ( ) 3 be the average tensor trace of ( ) and define the asphericity: Finally, the shape parameter reads: 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 50%. 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, 1 and 2 , are computed in decibels following the convention 10log 10 � 1 2 �. 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.

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: Then, the gelation ratio order parameter reads: Defining the largest cluster size as: we let the normalized largest cluster size order parameter be: The average cluster size reads: and the normalized cluster size order parameter is given by: The standard deviation cluster size is defined as: therefore, the normalized standard deviation cluster size order parameter is: The reasoning behind dividing gelation-related quantities by filaments is to map them on the segment [0,1].

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 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][62][63][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][66][67][68]. 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 4, 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: where |E| is the number of edges between node and node : ( , ) ∈ . is a symmetric adjacency matrix of : is the degree of node and equals the number of edges connected to . , is Kronecker delta and equals 1 if = and 0 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: where is the clustering coefficient of node : 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: 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: The eigenvector centrality of node is the -th element of the eigenvector defined by: where max ≥ 1 ≥ ⋯ ≥ |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: 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: 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: 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: 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: 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: 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: 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: We used the python KernelDensity estimator the Scikit Python library [75] to evaluate the selfinformation 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:

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 1 μm 3 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 ≥ 3 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.

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)(4)(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 (10dB) 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 (20dB) greater than those of low motor-filament ratio. , 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.

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 nonlocal 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 0.14, the average node degree is 140.7, and the average path length 6.1 whereas for the graph representing the hexavalent multilinkers ( = 6), the graph density of 0.42, the average node degree of 422.2, and the average path length of 3.1. 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= 1 s, 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)  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 ≥ 4 , (Fig. 9(b)) at valencies ≥ 4 , 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 0.7 to 0.5 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 0.14, an average node degree of 140.7, and an average path length of 6.1 compared with graph density of 0.42, an average node degree of 422.2, and an average path length of 3.1 in the system with hexavalent multilinkers.
At 1:100 motor-filament concentration ratio together with multilinkers with valency ≥ 3, 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 ≥ 4, high motor concentration enhances the assortativity reduction and causes arborization of actin bundles; motors reduce the assortativity of the network from 0.7 to 0.5 for the cases of hexavalent ( = 6) and heptavalent ( = 7) multilinkers. 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 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). 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 ≥ 4 as also shown for the 1:1 motor-filament concentration ratio in Fig. 8(d).

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 coarsegrained 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 actinbinding 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 ≥ 4 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