Topology determines force distributions in one-dimensional random spring networks

Networks of elastic fibers are ubiquitous in biological systems and often provide mechanical stability to cells and tissues. Fiber reinforced materials are also common in technology. An important characteristic of such materials is their resistance to failure under load. Rupture occurs when fibers break under excessive force and when that failure propagates. Therefore it is crucial to understand force distributions. Force distributions within such networks are typically highly inhomogeneous and are not well understood. Here we construct a simple one-dimensional model system with periodic boundary conditions by randomly placing linear springs on a circle. We consider ensembles of such networks that consist of $N$ nodes and have an average degree of connectivity $z$, but vary in topology. Using a graph-theoretical approach that accounts for the full topology of each network in the ensemble, we show that, surprisingly, the force distributions can be fully characterized in terms of the parameters $(N,z)$. Despite the universal properties of such $(N,z)$-ensembles, our analysis further reveals that a classical mean-field approach fails to capture force distributions correctly. We demonstrate that network topology is a crucial determinant of force distributions in elastic spring networks.


I. INTRODUCTION
Networks-and their topologies-have been studied in a broad range of disciplines, leading to terms like social, economic, biological, or chemical networks, and, of course, mechanical networks [1]. Here we focus on the latter and expand the theoretical and numerical analysis introduced in a companion short paper [2]. Networks of filamentous proteins, polysaccharides or nucleic acids, essentially all semiflexible filaments, play important roles for the mechanics and stability of biological cells and tissues [3,4]. An important design feature of biological materials is the response to large loads, including failure, rupture, damage limitation and their recovery properties. To understand failure that starts with the rupture of single filaments when the local force exceeds a threshold, it is crucial to understand force distributions in filament networks. It turns out that topology plays a critical role for the distribution of forces in elastic (e.g, polymer) networks, but this topic has received little attention to date.
The quantitative analysis of force distributions within random polymer networks has largely relied on computational modeling [5,6]. Analytical descriptions of fiber networks have primarily used effective-medium [7][8][9] or mean-field [6,10,11] approaches. Effective-medium 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 * These two authors contributed equally. † cfs@physik3.gwdg.de ‡ wardetzky@math.uni-goettingen.de 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 the effect of external force applied to the network, we employ a generation procedure that inserts springs with pre-strain so that the resulting initial configurations 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 arXiv:1707.01549v2 [cond-mat.soft] 7 Jul 2017  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 that are each uniformly distributed as U(−0.5, 0.5), but, since they are coupled by integer winding numbers, not mutually independent [12] as random variables.
We seek to characterize the length (i.e. force) distributions of springs in networks after they have relaxed to mechanical equilibrium. Relaxation preserves network topology, i.e., it preserves its graph together with a set of winding numbers, that arise 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 uniquely determines network topology and results in a known linear solution operator for the respective mechanical equilibrium. However, a network ensemble, with a given connectivity and number of nodes includes many topologies. This leads to a random solution operator, which makes it more difficult to determine the ensemble-averaged distribution of relaxed lengths.
Motivated by experiments, where explicit information on particular realizations is hard to obtain, we study ensembles with a fixed number of nodes N and average degree z, henceforth called (N, z)-ensembles. Surprisingly, such ensembles have well defined force distributions despite varying topologies. Explicitly accounting for these unknown underlying topologies makes our approach different from a mean-field description.

II. ANALYTICAL THEORY
Formally, as already described in [2], our model can be described 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 lengthsl and the signed cycle matrix C ∈ Z m×N z/2 , described below. 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. A fundamental cycle is defined as a cycle that occurs when adding a single edge to a spanning tree of the graph. There are N −1 edges in the spanning tree, so N z/2−(N −1) edges can be added. Therefore, there are N (z/2 − 1) + 1 fundamental cycles. Note that the choice of fundamental cycles corresponds to the choice of a basis and is therefore not unique. The solution to Eq. (1), however, is independent of this choice (Appendix A).
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.
It is noteworthy to point out that the problem presented above is equivalent to the classical problem of determining the currents (here l) in an electrical network. Force balance (or minimizing the energy l T l) is equivalent to Kirchhoff's current law (signed currents add up to zero at a node) and-assuming unit resistances-the cycle constraints (Cl = g) correspond to Kirchhoff's voltage law (voltages in a closed loop sum up to zero), where the winding numbers g j represent voltage sources.
Equation (1) defines a quadratic programming problem with a unique analytic solution: which can be explicitly computed for each realization via, e.g., the optimization library IPOPT [13].
To express the resulting force distributions of an (N, z)-ensemble we consider the expected histogram of the vector l * of random variables. This results in a univariate probability density for the final spring lengths. For a particular realization, the corresponding cumulative histogram H l * is given via where 1 A is the indicator function (one if A is true, zero otherwise). The quantity H l * ( * ) measures the number of elements in l * with values less than or equal to * . We are interested in a univariate cumulative distribution function (cdf) F l * , which we define as the expected value of the cumulative histogram of an (N, z)-ensemble: The quantity F l * ( * ) is the average over the marginal distribution functions of the individual l * i . This result defines the corresponding univariate probability density (expected histogram), i.e., By decomposing the final length vector l * into initial lengthsl and length changes ∆l, i.e., l * =l + ∆l, we compute and therefore with Eq. (5): Remember that the initial spring lengthsl i are identically distributed, i.e., pl i = pl. Equation (6) thus simplifies to: with p ∆l|l=¯ (∆ ) : Note that p ∆l|l=¯ , with the apparent dimensionality mismatch, is a shorthand notation that does not mean that l i =¯ for all indices i, but instead, corresponds to the average over all possible events thatl i =¯ for some index i. In this sense, the nth raw moment of the conditional probability density Eq. (8) is defined as follows: In the following we characterize the conditional probability density given in Eq. (8) that completely determines the final distribution of spring lengths given the initial distribution (Eq. (7)). Reconsidering Eq. (2), we write Equation (10) relates ∆l tol and a random matrix S, both of which vary with the topology of each realization. It is therefore challenging to obtain p ∆l|l=¯ explicitly, especially since the individuall i are not mutually independent. Instead, we consider the first two moments of the probability distribution, E(∆l|l =¯ ) and Var(∆l|l =¯ ), and investigate under which conditions ∆l|l =¯ is approximately normally distributed. In the following we will work with conditional random variables, so we now highlight two important aspects of our generation procedure that will be used extensively. The first is that each graph cycle, and therefore constraint, contains at least three edges, implying that the edge lengths are pairwise independent as random variables. The second aspect is that we can fix the abstract graph structure in our generation procedure, leading to (N, z)-ensembles with varying winding numbers, but with a constant S-matrix (e.g., Fig. 1). These fixedgraph-ensembles still contain identically, uniformly distributed random variablesl j ∼ U(−0.5, 0.5).

A. Conditional mean
In this section we compute E(∆l|l =¯ ) for (N, z)ensembles. We first derive the conditional mean for a fixed-graph-ensemble, i.e., E[(∆l|l =¯ )|S], and then generalize the result to (N, z)-ensembles. Equations (9) and (10) lead to: where we used the fact that fixed-graph-ensembles have uniformly distributed edge random variables that are pairwise independent. We further make use of our knowledge about the graph's cycle matrix C to determine tr S. First note that by definition tr S = tr P − N z/2 (Eq. (10)). The projector property of P (i.e., P 2 = P) leads to tr P = dim(Im P) = N z/2 − dim(ker P), because P has eigenvalues 0 and 1 only. Furthermore, ker P = ker C, by definition (Eq. (2)), and hence tr S = − dim(ker C). Recall that C contains N (z/2 − 1) + 1 linearly independent rows corresponding to a set of fundamental cycles of the graph, i.e., C has full rank and so dim(ker is an invariant of the (N, z)-ensemble as it surprisingly only depends on the number of nodes in the graph. Making use of this invariance together with the general property of the expected value, E(X) = E Y [E(X|Y )], we combine Eqs. (11) and (12) to obtain:

B. Conditional variance
The conditional variance Var(∆l|l =¯ ) remains challenging to express analytically for arbitrary z and N . To compute the conditional mean, we used the essential fact that expectation is always additive regardless of the dependencies between the random variables. This is not true for variances. The variance is only additive if the terms are pairwise independent. While the edge random variables are pairwise independent, they are in general not conditionally pairwise independent since the remaining two edge random variables of a triangle (cycle with three edges) with one edge length fixed are coupled by the fact that they sum to an integer winding number. For an (N, z)-ensemble, we do not know the abstract graphs, let alone their triangle structures, making the general computation of the conditional variance difficult.
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 with a known triangle structure. Both are symmetric (i.e., vertex-and edge-transitive [14]). In particular, edge-transitivity (informally: edges are indistinguishable from each other) allows us to reduce to a single entry in l * , since p l * = p l * i . The single component l * i is given by a weighted sum of identically distributed, but dependent random variables (Eq. (2)), which we analyze to derive Var(∆l|l =¯ ) explicitly. We first present the conditional variance derivation for these extreme cases and then discuss the more general intermediate-connectivity regime 2 < z < N − 1 that contains ensembles of multiple graphs. The complexity in this regime is highlighted by the intricacies involved in deriving the variance for a fixed graph (e.g., the complete graph), which already requires a special choice of basis to obtain a tractable expression for (CC T ) −1 .

Cycle graph
For the cycle graph (z = 2), there is only one cycle that contains all N edges. Therefore, the cycle matrix can be written as C = (1, 1, . . . , 1) ∈ R N . It follows that Eq. (2) simplifies to where g = N j=1l j is the winding number of the cycle and I ∈ R N is the vector of ones. We derive the conditional variance Var(∆l|l =¯ ) = N −1 N i=1 Var(∆l i |l i =¯ ) for the cycle graph. By edge-transitivity Var(∆l|l =¯ ) = Var(∆l i |l i =¯ ) and by Eq. (14), For the cycle graph, if N > 3, the conditional edge random variables are pairwise independent, so we compute: a constant that is independent of the initial spring length , showing that Var(∆l|l =¯ ) = El[Var(∆l|l)]. Conditional pairwise independence only holds for N > 3 because, for N = 3, if we condition on one length (l 1 =¯ ) the remaining two lengths are dependent viā l 3 = g −¯ −l 2 . Therefore, a similar computation for the conditional variance of a general ensemble does not hold, since each graph may contain triangles with conditionally pairwise dependent edges.

Complete graph
For the case of the complete graph (z = N − 1), 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 (i.e., spanning tree). This choice of basis leads to a tractable expression for (CC T ) −1 , which can then be applied to reformulate the problem in terms of conditionally independent winding number random variables.
We choose a spanning tree as shown in Fig. 3. For the following derivation, we label the edges such that the first N − 1 edges correspond to the edges of the spanning tree. The other m = N (N − 1)/2 − (N − 1) edges are the ones that are added to the spanning tree to construct the fundamental cycles (here: all triangles). We order the cycles in C according to these edges and decompose the cycle matrix into two parts: where I is the identity matrix. Our first result is that for spanning tree edges l * st := where (C T C) st corresponds to the first N − 1 rows of C T C. The importance of this result is that the symmetries of the complete graph allow it to extend to all edges. Indeed, each vertex of the complete graph defines a spanning tree as shown in Fig. 3; therefore every edge can be seen as such a spanning tree edge. Independence of the solution of the choice of the cycle basis, and therefore spanning tree, implies that the following derivations hold for all edges in the graph. Since, by Eqs. (2) and (17), proving Eq. (18) is equivalent to showing that We can construct A T A ∈ R (N −1)×(N −1) explicitly: The ijth entry counts the number of cycles that are shared by the edges i and j-with a contribution of 1 if the edges have the same orientation with respect to the cycle, and −1 otherwise. As can be seen in Fig. 3, two spanning tree edges only share one cycle with opposite orientations, where J ∈ R (N −1)×(N −1) is the matrix with ones everywhere. Substitution into Eq. (19) yields: Equation (21) holds true since each column of A T contains two nonzero entries, 1 and −1, due to all fundamental cycles (triangles) involving two spanning tree edges with opposite orientations (Fig. 3).
We have therefore shown that for the spanning tree edges of the complete graph, the following relation holds: where (C T g) st is the vector of the first N − 1 entries of C T g, and g = Cl is the vector of winding numbers (Eq. (1)). Equation (22) Observe that the g j |l st=¯ are independent random variables since their only potential dependence, their common edge, is conditioned out. Each winding number g j |l st=¯ corresponds to a fundamental cycle that is a triangle (Fig. 3), i.e., involves only three edges of which one is fixed. Therefore we cannot use conditional pairwise independence of edge lengths as in the case of the cycle graph to compute the variance. Instead, we derive the winding number distribution of a triangle explicitly.
The initial edge lengths are distributed asl j ∼ U(−0.5, 0.5), so the winding numbers can only attain three values {−1, 0, 1}. In particular, g j |l st=¯ =¯ +l j1 + l j2 , where we choose positive signs since the lengths are distributed symmetrically around zero. We compute the probability of g j |l st=¯ attaining the value zero: where χ [−0.5,0.5] (·) is the characteristic function on the interval [−0.5, 0.5]. The remaining probability is assigned to either g j |l st=¯ = 1 or g j |l st=¯ = −1 depending on whether the given¯ is positive or negative: Using Eq. (23), conditional independence of the g j |l st=¯ , and their probability distribution, Eqs. (24) and (25), we have: In contrast to the cycle graph (Eq. (16)), the conditional variance Var(∆l|l =¯ ) = Var(∆l st |l st =¯ ) for the complete graph (Eq. (27)) depends on the initial spring length , as is shown in Fig. 2.

Intermediate-connectivity regime
For the intermediate-connectivity regime, 2 < z < N −1, a tractable expression for (CC T ) −1 , as for the complete graph, remains elusive; however, numerical data suggest that the 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 [12] states: We first compute Var(∆l), again, by initially fixing S and considering Var(∆l|S). With Eqs. (9) and (10) we compute where the second equality follows from fixed-graph ensembles having uniformly distributed edge random variables that are pairwise independent. Again, we use that P 2 = P, hence S 2 = −S, and therefore N z/2 j=1 S 2 ij = −S ii . Insertion into the equation above yields: where the second equality is due to Eq. (12). Application of the law of total variance gives: where the second term in the sum vanishes using Eq. (10) since E(l i ) = 0 (analogous to the computation in Eq. (11)). From Eq. (13) we can use Varl[E(∆l |l)] = (2/z(1 − 1/N )) 2 Var(l) and therefore find by substituting into Eq. (28): If ∆l|l =¯ were normally distributed, having estimates for mean and variance (Eqs. (13) and (31)) 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 . 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 [15], where independence holds beyond a certain time window, in our case 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 [16,17] (a deviation-bound version of the central-limit theorem) to obtain a quantitative bound on the distance to a normal distribution.
The rest of this section is split into three parts. We begin by proving the results for the cycle and complete graph, and then investigate the intermediateconnectivity regime. Throughout, note how the intricacies of the proofs of the extreme cases are further complicated in the intermediate-connectivity regime, where ensembles have varying graph structure and lack symmetry.

Cycle graph
For the cycle graph (z = 2), we prove that ∆l i |l i=¯ is normally distributed in the limit N → ∞. The key idea is a relaxation of the integer constraint to an interval constraint. Using ∆l i = g/N −l i = N j=1l j /N −l i (Eq. (14)) we introduce the standardized (E(Y N ) = 0, Var(Y N ) = 1) random variable and compare its cdf F Y N (x) to that of the standard normal Φ 0,1 (x). We find that where we have used the fact that g ∈ Z and · denotes the floor operator. Using g| li=¯ =¯ + N j=1, j =il j we can formalize the integer relaxation by expressing the probability of the conditional winding number as follows: where ν N −2 (t) corresponds to the probability density of the sum of N − 2 uniformly distributed independent random variables on the interval [−0.5, 0.5]. We call this random variable U N −2 . There are only N − 2 independent random variables because one of the N lengths is fixed to¯ and another one is determined to make sure that an integer winding number is attained for g. Substituting Eq. (38) into Eq. (37) and expressing with the random variable δ(x) ∈ [0, 1), leads to: We are interested in the distance to the cdf Φ 0,1 (x) of the standard normal distribution. To calculate this distance, this we perform a change of variables which results in a standardized sum of uniforms U N −2 /(σl √ N − 2), where σl = Var (l). However, this causes the cdf of our random variable and the standard normal cdf to have different arguments. We therefore split the computation into two steps: one that measures the distance to a shifted standard normal cdf, and the other that measures the deviations introduced by this shift. With the shorthand notation ξ(x) := x (N − 1)Var(l) + 0.5 − δ(x), the described procedure corresponds to the following computation: .
I(x) can be bounded using the Berry-Esseen theorem.
Bounding II(x) requires a detailed case analysis (see Ap-pendix B for details). We arrive at: Therefore, the cdf of ∆l i |l i=¯ converges to a normal distribution with the rate (N − 2) −1/2 , independent of¯ . Since we showed that F ∆li|li=¯ is independent of the edge i, Eq. (4) implies F ∆l|l=¯ (x) = F ∆li|li=¯ (x), and therefore F ∆l|l=¯ converges to a normal distribution as well.

Complete graph
For the complete graph, our proof of normality relies on the reduction to spanning tree edges as outlined in the conditional variance section. In particular, this allows us to write the cdf in terms of winding number random variables that are all triangles that share a common edge. Conditioning on this edge then yields independence, not of the length variables, but of these winding number random variables, which allows us to apply the Berry-Esseen theorem.
To measure how far ∆l st |l st =¯ is from being normally distributed for finite N , we look at the standardized random variable and compare its cdf to the one of the standard normal. Using Eq. (23) and the probability distribution of g j |l st =¯ , Eqs. (24) and (25), we obtain: and therefore All {g j |l st=¯ } N −2 j=1 are independent since the corresponding cycles only share one edge, which is the one that we condition on. We can thus apply the Berry-Esseen theorem (Appendix B) to show that for N ≥ 3, with C < 0.4748, ρ = E( g j |l st=¯ −¯ 3 ), and σ 2 = Var(g j |l st=¯ −¯ ) = |¯ | −¯ 2 . Computing ρ = (|¯ | −¯ 2 )(¯ 2 + (1 − |¯ |) 2 ) via Eqs. (24) and (25) we arrive at for |¯ | > 0. This proves convergence of the cdf of ∆l i |l i=¯ to a normal distribution with the rate (N − 2) −1/2 . For¯ = 0, ∆l i |l i=¯ ∼ δ 0 (Dirac delta distribution around zero), it can only attain the value zero because P (g j |l st=¯ = 0) = 1. Note the¯ -dependence in Eq. (50), which is in stark contrast to the¯ -independent bound for the cycle graph (Eq. (41)). For the complete graph, the approximation with a normal distribution becomes worse as¯ approaches zero.

Intermediate-connectivity regime
Recall that in the intermediate-connectivity regime, 2 < z < N − 1, the (N, z)-ensembles contain graphs with varying 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. 4).

D. Density approximation
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. (13) and (31). Using Eqs. (7) and (51), we obtain an explicit representation for the final length distribution p l * ( * ) in mechanical equilibrium (Appendix C). In Fig. 5 we compare this analytical expression to ensembles of simulated networks and observe excellent agreement.

III. COMPARISON TO A MEAN-FIELD APPROACH
In order to evaluate the significance of our graphtheoretical analysis we compare it to a mean-field (mf) approach which neglects all topological features other than the local degree of connectivity. In contrast to the graph-theoretical model, where z refers to the average degree of a node, the mean-field approach assumes that each node is connected to exactly z other nodes. Moreover, the node displacement u node during relaxation is calculated as if all other nodes in the network were fixed. Therefore u node = f node /z, where f node = z i=1l i is the initial force acting on the node via the springs attached to it. The displacement ∆l of a spring is given by the difference of the displacements of the two nodes that arel FIG. 6. Mean-field approach (here: z = 4): The edgel connects two nodes (n1, n2). The initial node forces (gray arrows) are given by fn 1 = z−1 i=1l n 1 ,i +l and fn 2 = z−1 i=1l n 2 ,i −l. The springl contributes to the forces with different signs because its length is measured from n1 to n2 (depicted by the black triangle). While displacing an individual node by un k = z −1 fn k introduces force balance at that node, this approach neglects that nodes/edges are coupled, i.e., force balance has to be established at all nodes simultaneously, as in Eq. (1).
connected by this edge: where we have taken into account that the two nodes share one spring, namelyl (Fig. 6). All springs are assumed to be independent identically distributed random variables with mean zero and variance Var(l). For the conditional mean, we have with Eq. (52): The mean-field result agrees with the exact solution Eq. (13) in the limit N → ∞, i.e., there is no significant difference for large node numbers. In contrast, we will show at the end of this section that for the variance, the mean-field solution differs substantially from the exact result, even in the limit N → ∞. When considering normality of ∆l|l =¯ , the mean-field approach allows us to directly apply the Berry-Esseen theorem (Appendix B) because all edges are treated as independent. Defining the normalized random variable the theorem implies: The mean-field approach yields convergence to a normal distribution with a rate proportional to (z − for which we showed convergence to a normal distribution even though z is constant. In the intermediateconnectivity regime, both the mean-field as well as our graph-theoretical approach suggest that ∆l|l =¯ can be approximated by a normal distribution. To complete the evaluation of our approach, it is therefore critical to also compare the second moments, i.e., the variances of the mean-field and graph-theoretical approach. For the unconditional variance, we obtain using Eq. (52): Var(∆l)| mf = z −2 4Var(l) + 2(z − 1)Var(l) (57) Clearly, this expression does not agree with the exact graph-theoretical value Eq. (30), even in the limit N → ∞. A mean-field approach assumes the conditional variance is constant and therefore equal to its expected value El[Var(∆l|l)]| mf = Var(∆l|l =¯ )| mf = 2/z(1 − 1/z) Var(l). For the cycle graph, we proved the conditional variance is indeed constant (Eq. (16)). However, we showed that the other extreme, the complete graph, exhibits non-constant conditional variance (Eq. (27)). For (N, z)-ensembles in the intermediateconnectivity regime, we observe a continuous transition between the two extremes (Fig. 7). Therefore, for the biological regime (z 4), we approximated the conditional variance with its constant expected value El[Var(∆l|l)] (Eq. (31)). However, it is exactly the regime z 4 where the graph-theoretically derived expected conditional variance El[Var(∆l|l)] and the mean-field quantity El[Var(∆l|l)]| mf exhibit the largest discrepancy (Fig. 7).

IV. DISCUSSION AND CONCLUSIONS
In conclusion, we have presented a probabilistic theory of force distributions in one-dimensional random spring networks on a circle. Here we have 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 length changes conditioned on an initial configuration. For the two extreme cases, the cycle graph and the complete graph, we could prove convergence of this distribution to a normal distribution. A systematic analytical treatment of the-less symmetricintermediate-connectivity regime 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 present 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 correctly reproduce the mean and the variance of the relevant distributions. The error is particularly pronounced for the-biologically 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 of network topology 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.

ACKNOWLEDGMENTS
The authors would like to thank Friedrich Bös, Alexander Hartmann, and Fabian Telchow for fruitful discussions. Funding from the Deutsche Forschungsgemein-schaft (DFG) within the collaborative research center SFB 755, project A3, is gratefully acknowledged. C.F.S was additionally supported by a European Research Council Advanced Grant PF7 ERC-2013-AdG, Project 340528.
Appendix A: Independence of the choice of cycle basis A change of cycle basis corresponds to the transformatioñ where Q ∈ GL(m) is an arbitrary change of basis matrix for the cycle space, and U ∈ O(N z/2) is a permutation matrix that corresponds to relabeling the edges of the graph. The independence of the solution of the cycle matrix means that, given the change of basis in Eq. (A1) and the solution (Eq. (2)) to the transformed problem,l * = U l * .

Appendix B: Upper bounds for I(x) and II(x)
For a uniform upper bound on the first term I(x), we can apply the Berry-Esseen theorem, which is stated as follows [18].