Topology counts: force distributions in circular spring networks

Filamentous polymer networks govern the mechanical properties of many biological materials. Force distributions within these networks are typically highly inhomogeneous and, although the importance of force distributions for structural properties is well recognized, they are far from being understood quantitatively. Using a combination of probabilistic and graph-theoretical techniques we derive force distributions in a model system consisting of ensembles of random linear spring networks on a circle. We show that characteristic quantities, such as mean and variance of the force supported by individual springs, can be derived explicitly in terms of only two parameters: (i) average connectivity and (ii) number of nodes. Our analysis shows that a classical mean-field approach fails to capture these characteristic quantities correctly. In contrast, we demonstrate that network topology is a crucial determinant of force distributions in an elastic spring network.

Filamentous polymer networks are ubiquitous in nature.They make up the cytoskeleton of animal cells and form the scaffold of the extracellular matrix in, e.g., connective tissue.These networks determine the mechanical response of cells and tissues and support elastic forces under external or internal loading at both mesoscopic and macroscopic scales.The force distributions within such networks can be highly inhomogeneous [1,2].Internal forces in these typically non-equilibrium networks result mostly from molecular motors [3][4][5] in the cell cytoskeleton or, on a larger scale, from cells embedded in extracellular matrices, such as platelets in blood clots [6,7].
The quantitative analysis of force distributions within random polymer networks has largely relied on computational modeling [1,8].Analytical descriptions of filamentous networks have primarily used effective-medium [9][10][11] or mean-field [8,12,13] approaches.Effectivemedium theories rely on mapping a disordered system to an ordered one.It is unclear, however, how force distributions change under this mapping.Mean-field approaches do not consider the full network topology, but only the local degree of connectivity.We show that such an approach fails to describe force distributions even for a very simple model system; in fact, topological features, i.e., cycles/loops in the networks, cause global coupling that remains prevalent even when the system becomes large.
The simple model system that we consider here consists of ensembles of one-dimensional random spring networks on a circle.Considering such networks is equivalent to applying periodic boundary conditions in one dimension.To model a generically forced system we employ a generation procedure that results in initial configurations that are not in mechanical equilibrium.We then study the resulting force distributions of the relaxed systems.
We generate initial network configurations as follows (Fig. 1): (i) Place N node positions (indexed from 1 to N ) drawn from a uniform distribution on the circle.(ii) Connect these nodes in the order given by their indices into one connected cycle via springs.We always connect consecutive nodes via the shorter of the two possible distances.Note that the cycle may wrap around the circle zero, one, or multiple times (Fig. 1 (c)).This step guarantees that each network will always have only one connected component and prevents dangling ends.
(iii) Connect further node pairs randomly, such that each node pair is connected by at most one spring, until the network contains N z/2 springs, where the average degree of connectivity z is chosen such that N z/2 is an integer.Each spring is linear, has rest length zero, and unit spring constant.Its length is measured along the circumference of the circle.In order to encode this construction in an unambiguous manner we work with signed spring lengths as degrees of freedom.The orientation of a spring is chosen such that it goes from a node of lower index to a node of higher index.This is an arbitrary choice, but defined orientations are essential in our formalism.The sign of the spring length is chosen to be positive if its orientation on the circle points counter-clockwise and negative otherwise.
The network can be encoded within a graph representation, where the springs together with their orientations are the directed edges of the graph, with signed lengths as edge weights (Fig. 1 (b)).To lie on the circle, the graph and edge weights must be compatible in the sense that the sum of the edge weights around each cycle of the graph is equal to an integer, which we refer to as its winding number g.Our network generation procedure guarantees this compatibility.It results in a random directed Hamiltonian graph, i.e., a graph that contains a cycle that visits each node exactly once, with N nodes and average degree z.This graph comes equipped with compatible initial spring lengths/edge weights { li } The network contains two fundamental cycles, for example: {l1, l2, l3, l4} and {l4, l5, l6}.After choosing arbitrary orientations for both cycles (gray arrows), we construct linear constraints that fix their winding numbers (Eq.( 1))here: l1 + l2 + l3 − l4 = −1 (winds around circle once) and l4 + l5 + l6 = 0 (contractible).(c) The abstract cycle graph (z = 2) with N = 5 (left) and three realizations on the circle with distinct topologies (same graphs but different winding numbers g).Top and bottom row show initial and corresponding relaxed configurations, respectively.Note that, for visualization purposes, overlapping springs are drawn with a slight offset.
since they are coupled by integer winding numbers, not mutually independent [14] as random variables.
We seek to characterize the length (i.e.force) distributions of springs in networks after they have relaxed into mechanical equilibrium.Relaxation preserves a network's topology, i.e., its graph together with a set of winding numbers, that arises from the generation process.Note that networks sharing the same graph may have different sets of winding numbers, and therefore distinct relaxed states (Fig. 1 (c)).A particular realization of an initial network (see above) uniquely determines network topology and results in a known linear solution operator for the respective mechanical equilibrium.However, a network ensemble realizes many topologiesyielding a random solution operator-making it more difficult to determine the ensemble-averaged distribution of relaxed lengths.
Motivated by experiments, where explicit information on particular realizations is hard to measure, we study ensembles with a fixed number of nodes N and average degree z. Surprisingly, such ensembles have characterizable force distributions despite varying topologies.Explicitly accounting for these unknown underlying topologies makes our approach different from a mean-field description.
Formally, our setup can be written as the following optimization problem: where l ∈ R N z/2 is the vector of all spring lengths and g ∈ Z m is the vector of winding numbers, which is determined by the vector of initial spring lengths l and the signed cycle matrix C ∈ Z m×N z/2 , described below.Note that all the above quantities are random variables.The first part in Eq. ( 1) minimizes the total elastic energy of the system, whereas the second part preserves the topology of the network by fixing the winding numbers of a set of m = N (z/2 − 1) + 1 fundamental cycles [15].Note that the choice of fundamental cycles corresponds to a choice of basis and is therefore not unique.The solution to Eq. ( 1), however, is independent of this choice [16].
After choosing a cycle basis, the C-matrix is constructed by specifying an orientation for each fundamental cycle and then setting C ji equal to: 1 if spring i is part of the jth fundamental cycle and their orientations agree, or −1 if their orientations are opposite, and 0 otherwise.For the example in Fig. 1 (a), the cycle matrix and vector of winding numbers are given by C 1 = (1, 1, 1, −1, 0, 0), C 2 = (0, 0, 0, 1, 1, 1), and g = (−1, 0) T , respectively.Note that winding numbers correspond to the signed number of times a cycle wraps around the circle.Contractible cycles have winding number zero.If all cycles were contractible, then Eq. ( 1) would have a trivial solution with all springs collapsed to a single point; it is only the presence of nontrivial cycle constraints that prevents this outcome.
Equation ( 1) defines a quadratic programming problem with a unique analytic solution.Written in terms of the spring length changes ∆l during relaxation to the final configuration l * , the solution is: which can be explicitly computed for each realization via, e.g., the optimization library IPOPT [17].
To express the resulting force distributions of our ensembles we consider the expected histogram of the vector l * of random variables.This results in a univariate probability density p l * , which is given by the average of the individual spring densities, i.e., p l * ( * ) := 2 N z N z/2 i=1 p l * i ( * ) (for details see [16]).Using l * i = li +∆l i (Eq.( 2)) we compute for each component: Since the initial spring lengths li are identically distributed, i.e., pl i = pl, we obtain that with In the following we characterize the conditional probability density Eq. ( 5) that completely determines the final distribution of spring lengths given the initial distribution (Eq.( 4)).Equation ( 2) relates ∆l to l and a random matrix S, which vary with the topology of each realization.It is therefore challenging to obtain p ∆l| l= ¯ explicitly, especially since the individual li are not mutually independent.Instead, we consider the first two moments, E(∆l| l = ¯ ) and Var(∆l| l = ¯ ), and investigate under which conditions ∆l|l = ¯ is approximately normally distributed.
Equations ( 2) and ( 5) lead to (see [16] for the derivation): where tr S = 1 − N is an invariant of the ensemble that surprisingly only depends on the number of nodes in the graphs, not on their respective topologies.We compare Eq. ( 6) to a mean-field (mf) approach, where each node is displaced as if all other nodes in the network were fixed.In this case E(∆l| l = ¯ )| mf = −2 ¯ /z [16]; in particular, the mean-field result agrees with the exact solution Eq. ( 6) in the limit N → ∞, i.e., there is no significant difference for large node numbers.In contrast, we will show that for the variance, the mean-field solution differs substantially from the exact result, even in the limit N → ∞.
The conditional variance Var(∆l| l = ¯ ) remains challenging to express analytically.Indeed, in general, there are many graphs realizing the same z and N , each with its own topology that may introduce nonzero covariance between the edge lengths.
For two extreme cases, namely the cycle graph (z = 2, N > 3) and the complete graph (z = N − 1, each node connected to every other node), there exists only a single possible graph, respectively, each being symmetric (i.e., vertex-and edge-transitive [18]), allowing us to derive Var(∆l| l = ¯ ) explicitly.In particular, edge transitivity allows us to reduce to a single entry in the l * vector, which is given by a weighted sum of identically distributed, but dependent random variables (Eq.( 2)).
For the case of the cycle graph, an entry in Eq. ( 2) simplifies to ∆l i = N −1 j =i lj .Using conditional pairwise independence of the initial edge random variables allows for direct computation of Var(∆l| l = ¯ ) = (N − 1)/N 2 Var( l).
For the case of the complete graph, the derivation of the conditional variance is significantly more involved.In order to obtain manageable algebraic expressions, one needs to carefully choose the cycle basis.This choice is detailed in [16] and leads to a tractable analysis of (CC T ) −1 , which can then be applied to reformulate the problem in terms of conditionally independent winding number random variables, leading to Var(∆l| For the intermediate-connectivity regime, 2 < z < N − 1, a similar approach remains elusive; however, numerical data suggest that the conditional variance exhibits a continuous transition between the two extremes (Fig. 2).We also observe that the conditional variance is approximately constant given that z N .This is the most relevant case for biological networks where typically z 4. For z N , we may thus approximate Var(∆l| l = ¯ ) ≈ El[Var(∆l| l)], which we now derive.
The law of total variance [14] states: From Eqs. ( 2) and ( 5) it follows that where we have used that S 2 = −S (see [16] for details).
Clearly, the expressions do not agree in the limit N → ∞.In particular, for sparsely connected networks (small values of z), there are significant deviations, independent of the number of nodes in the network.Using Eq. ( 6) we also have that Varl and therefore by substituting into Eq.( 7): Now, if ∆l|l = ¯ were normally distributed, having estimates for mean and variance would be sufficient to fully characterize p ∆l| l= ¯ .Indeed, for the two extremes, cycle and complete graph, we can prove that ∆l|l = ¯ is normally distributed in the limit N → ∞, with a rate of convergence proportional to (N − 2) −1/2 [16].
This result might look like a direct application of the classical central limit theorem.However, since the edge lengths are not independent as random variables, more sophisticated techniques are required to represent the solution in terms of a suitable set of mutually independent random variables.In contrast to situations in time series analysis [19], where independence holds beyond a certain time window, the cycle constraints prohibit localization of dependencies.To deal with this problem, we reduce the number of variables by relaxing each integer cycle constraint to an interval constraint.Harnessing the resulting independence then requires a non-standard transformation of random variables, which complicates a direct application of the Berry-Esseen theorem [20,21] (a deviation-bound version of the central-limit theorem) to obtain a quantitative bound on the distance to a normal distribution.
Recall that in the intermediate-connectivity regime, 2 < z < N −1, the ensembles contain graphs with varying FIG. 4. Probability density p l * ( * ) for the final spring lengths for networks with N = 1000 and varying z.Solid black lines show the analytic expression for p l * ( * ) [16]; data points correspond to averages over 50 simulations.The error bars correspond to the standard deviation.For comparison, we show the initial uniform spring length distribution pl( ¯ ) as a gray dashed line.
cycle structures making a similar analysis significantly more challenging.In simulations, however, we observe that ∆l|l = ¯ is approximately normally distributed if z is sufficiently large (Fig. 3).
Our empirical observations and theoretical discussion above justify the following approximation for 3 ≤ z N : with the expressions for E(∆l| l = ¯ ) and El[Var(∆l| l)] given in Eqs. ( 6) and ( 9).Using Eqs. ( 4) and (10) we obtain an explicit representation for the final length distribution p l * ( * ) in mechanical equilibrium (see [16]).In Fig. 4 we compare this analytical expression to ensembles of simulated networks; we observe excellent agreement.
In conclusion, we presented a probabilistic theory of force distributions in one-dimensional random spring net-works on a circle.Here we regarded networks with initially unbalanced forces that relax into mechanical equilibrium.When drawing the analogy to a biological network, our approach, which focuses on the relaxation of the system after non-equilibrium starting conditions, is equivalent to assuming a separation of time scales where internal or external non-equilibrium processes slowly create forces in the network that rapidly equilibrate.
We developed a graph-theoretical approach that allows us to exactly compute mean and expected variance of the distribution of conditional length changes that completely determines the final length distribution.For the two extreme cases, the cycle graph and the complete graph, we could additionally prove convergence of this distribution to a normal distribution.A systematic analytical treatment of the-less symmetric-intermediate regime of connectivity is more demanding and not provided here.However, our results suggest an approximation that shows excellent agreement with simulation for the biologically relevant regime of connectivity, 3 ≤ z N .It is straightforward to generalize the approach we presented here to higher spatial dimensions d if the probability densities pl k for the components of the initial spring vectors are independent.In that case, due to the linearity of spring forces with extension, the optimization problem decouples into the spatial components.The probability density for the final spring vectors then is simply given as the product of the one-dimensional results: Hence, our results carry over to two-and threedimensional networks, which are more commonly studied in practice and are of biological and physiological relevance.Interestingly, a classical mean-field approach fails to capture the mean and the variance of the relevant distributions.The error is particularly pronounced for thebiologically most relevant-regime of low degrees of connectivity, and does not vanish in the limit of infinite node number.Our work demonstrates that network topology-here manifested as cycle constraints-is crucial for the correct determination of force distributions in an elastic spring network.
This opens the door for future research on the role that network topology plays in more complex elastic networks, e.g., in the presence of dynamics, spring nonlinearities or rupture.Moreover, the mixture of probabilistic and graph-theoretical techniques may prove useful for other types of network theories.

FIG. 2 .
FIG.2.Normalized conditional variance Var(∆l| l = ¯ )/El[Var(∆l| l)] as a function of ¯ for graphs with N = 100 and varying z.For each value of z, data points correspond to ensemble averages (repeated simulations) with 4.95 × 10 6 springs in total.We use local linear regression with 3 × 10 4 nearest neighbors to estimate the variance for different values of ¯ .The solid lines correspond to the analytically derived expressions for the cycle and the complete graph (illustrated in the insets).In the intermediate regime of connectivity, the variance shows a continuous transition between the two extreme cases.