Computational study of the geometric properties governing the linear mechanical behavior of fiber networks

Materials whose microstructure is formed by random fiber networks play an important role both in biology and engineering. So far, it still remains unclear which geometric properties of the fiber network determine the macroscopic mechanical properties of such materials. This paper presents a computational study based on a large number of representative volume elements of random fiber networks. Our study reveals that the linear mechanical properties of fiber networks (i.e., Young ’ s modulus and Poisson ’ s ratio) are largely determined by only four scalar key descriptors. These are the number of fibers per volume, the mean node valency, the mean fiber length, and the mean direction cosine between fibers adjacent to the same node. Number of fibers per volume and node valency were found to be responsible for around 80% of the variance of the mechanical properties, making them the two by far dominant microstructural descriptors. In the part of the configuration space covered by our study, we observed a linear or quadratic relationship between the above four scalar microstructural descriptors and the Young ’ s modulus. For the number of fibers per unit volume we propose a theoretical explanation for this simple


Introduction
Fibrous materials are abundant both in nature (e.g., soft tissues such as ligaments and skin or wood), and in engineering (e.g., textile materials, paper). Therefore, the development of numerical models to study their mechanical properties under divers mechanical loading conditions has received significant attention, in materials science, mechanical engineering, biomechanics and biophysics [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18]. The key to computational modeling of fibrous materials is often the ability to set up a socalled representative volume element (RVE) of their microstructure. The geometric features of these RVEs need to agree with the ones of the microstructure of the real materials with respect to relevant statistical descriptors [19][20][21][22][23]. If so, computational studies with the RVE can help understand the mechanical properties of the materials on the macroscale and also how certain changes of their microstructure could be used to optimize these.
The microstructure of fibrous materials is often defined by a random network structure. Over the last decades, numerous stochastic algorithms were developed and applied to construct RVE of random heterogeneous media [24][25][26][27][28]. For detailed reviews of the current state of research in microstructure characterization, reconstruction approaches and multi-scale modeling of heterogeneous materials the reader is referred to [29][30][31][32]. Generally, reconstruction of heterogeneous materials is an optimization problem, which, due to its complexity, typically requires an iterative numerical solution. Among various numerical approaches, simulated annealing (SA) has received much attention. For SA, one defines a desired random microstructure by a set of so-called descriptors (each of which characterizes a specific geometric property of the microstructure). Then one typically starts with some random initial configuration and uses SA to adjust it by random steps until obtains a microstructures whose descriptors match the desired ones. SA was introduced for the reconstruction of dispersions of particles based on correlation functions [33]. Afterwards, the method was extended to the reconstruction of general random heterogeneous media [34], thereby becoming applicable to multidimensional, multiphase and anisotropic structures. Since SA is a method that can be adapted quite easily to a specific problem, many studies have been performed in the past on the basis of SA [35][36][37][38].
Despite the many computational studies that have been performed with fiber networks, it remains to date still unclear which descriptors of their microstructure exactly govern their mechanical properties on the macroscale. To address this question, we propose in this paper an algorithm for the construction of random fiber networks. We use this algorithm to construct a large variety of random RVEs and compute their mechanical properties by finite element simulations to identify the descriptors that mainly govern the mechanical properties of fiber networks on the macroscale.
The present article is organized as follows. In Section 2, we briefly outline a large variety of different descriptors that can be used to characterize the geometry of fiber networks. In Section 3, we describe our algorithm to construct tailor-made RVE whose descriptors match the desired ones. In Section 4, we present the results of a computational study that reveal which of the large variety of considered descriptors are geometrically and mechanically most relevant. Finally, we summarized our results in Section 5.

Descriptors
Fiber networks consist of fibers connected to each other at so-called nodes. In order to characterize them or compare them with each other, one can use so-called descriptors. Each descriptor characterizes a different aspect of the geometry of the network. As we study herein random fiber networks, the descriptors are in general related to statistical properties of the networks. The major goal of this paper is to understand which descriptors are primarily governing the mechanical properties of fiber networks. In general there is an infinite number of possible descriptors. As rigorous mathematical proofs about the relevance or irrelevance of a descriptor appear to date impossible in most cases, we have to adopt a heuristic strategy. That is, in the following we present a finite set of the most common descriptors for networks from Euclidean geometry and graph theory. The role of each of these descriptors will then be examined in a computational study. This way, we cannot make mathematically rigorous statements about the set of mechanically relevant descriptors. However, we can at least address this question from a heuristic perspective that can be hoped to be sufficient especially for most purposes in materials research.

Morphological descriptors
Morphological descriptors are based on the geometry and topology of the network in the Euclidean space. With respect to fiber networks, the following ones are most prominent.
Node density ρ node is the number of nodes (vertices) in the RVE N node divided by the volume of the RVE. Valency distribution p v is the probability distribution of valency across the nodes in the RVE. The valency of a node v is the number of fibers (edges) connected to that node. Fiber length distribution p l is the probability distribution of the Euclidean length of the fibers l in the RVE. Direction cosine distribution p c is the probability distribution for the cosine c of the angles between all pairs of fibers connected to the same node. It describes how much the orientation of fibers joining at the same node is correlated. Connectedness of the networks is a boolean descriptor which declares whether the fibers in the RVE are all connected to each other in some way, that is, whether between any pair of nodes in the network there exists at least one connection path along the fibers (edges) of the network. Pore-size distribution function p pore (r) describes in the fiber network the probability for a point in the void phase that the nearest point in the fiber phase is located at a distance of at least r (see also the more general discussion in [39]). Radial distribution function p r (r) is a function of distance which describes in a system of point-like particles the probability of finding particles in the distance of r from a reference particle [33]. In this work we apply this descriptor in two different ways to characterize the solid phase (fiber phase). First, we consider the nodes the relevant set of particle and compute the probability of finding nodes in the distance of r from each other denoted as p r− node (r). Second, p r− fiber (r) which is the probability of finding fiber segments in the distance of r from each other. To this end, we discretize all the fiber within the RVE into segments. Then we compute the distance of all pair of segments within the RVE based on their central points. Geometric moment invariants I are very common and powerful tools for the recognition of objects and patterns in image processing. For an image described by the scalar intensity function f(x, y, z) in three dimensions, the geometric moments are [40] where the exponents i, j and k are non-negative integers and r = i +j +k is the order of the moment. The low order moments have physical concepts. m 000 can be interpreted as the mass of an object, x = m 100 /m 000 , y = m 010 /m 000 and z = m 001 /m 000 respectively as the coordinates of the center of mass in x-, y-and z-direction, and the second order moments m 200 , m 020 , m 002 as the moments of inertia around the x-, y-and z-axis respectively. Accordingly, the so-called central geometric moments can be defined as [40] By the combination of different geometric moments, we can define various geometric moment invariants which are unchanged under a group of transformations such as translation, rotation and scaling [41]. A large list of geometric moment invariants can be found in [42]. The invariants of higher order are often extremely small compared to the invariants of lower order and are often neglected in image processing. Thus also we herein use only the first three invariants In our fiber networks, we assumed f to be a Dirac-type function with infinite intensity on the (infinitesimally thin) fibers and zero intensity everywhere else. For the practical calculation of the geometric moment invariants, the fibers were divided into a finite number of segments, each associated with a weight according to the segment length, and then an approximation of the moment invariants was calculated based on this discretization.

Graph descriptors
In mathematics, graphs are formed by a set of vertices (nodes) and edges (links or lines) connecting them. So-called simple graphs are graphs where no vertex is connected by an edge to itself and where the connection between two vertices is always established by exactly one (rather than in general several) edges. Undirected graphs are graphs where the connections between vertices are bidirectional. Apparently, one can interpret networks of thin fibers (on which we solely focus herein) from an abstract point of view as graphs. Doing so, the fibers play the role of the edges and the nodes the one of the vertices, and the mathematical quantities that are usually applied to characterize graphs (graph descriptors) can be used to characterize also fiber networks. There are many different descriptors in graph theory [43][44][45][46], and it is not possible to study all of them in this paper. Rather we focus on a set of the most common ones briefly summarized in the following.
Clique number is the number of vertices in the largest clique of a graph. A clique is a subset of a graph in which there is an edge between each pair of vertices [45]. Domination number of a graph is the number of vertices in the smallest dominating set of that graph. A dominating set of a graph is a subset of the graph such that every vertex outside this subset is a neighbor of at least one vertex within the subset [43]. Independence number is the number of vertices in the largest independent set of a graph. An independent set of a graph is a subset of vertices within which no pair of vertices is directly connected by an edge [45]. Chromatic number of a graph is the smallest number of colors needed to color the vertices of that graph such that the color of any two adjacent vertices is different [45]. Clustering coefficient in graph theory is a tool to indicate the tendency of the nodes to cluster together. It can be defined either for each node individually (local clustering coefficient) or for the whole graph (global clustering coefficient). The local clustering coefficient of a specific node is the ratio of the actual number of edges between its neighbors and the maximal possibly number of edges between its neighbors. The global clustering coefficient is the ratio of the number of closed triplets over the total number of triplets in the network. A triplet is a set of three vertices, at least one of which is connected to the other two. If also the other two are connected to each other, the triplet is called a closed triplet. The global clustering coefficient is well-known to be a measure of the clustering in the network [46] S-metric of a graph is the product of the valencies of all vertices [44].
Maximum eigenvalue of adjacency matrix: the adjacency matrix of a graph is a square matrix whose size equals the number of vertices. Its ij-elements are the numbers of edges connecting the i-th and j-th vertex. For the simple undirected graphs corresponding to fiber networks the adjacency matrix is a symmetric matrix of zeros and ones where the diagonal elements are all zero. Thus the adjacency matrix has real eigenvalues. Its maximum eigenvalue is a frequently used descriptor to characterize graphs. [43] Maximum eigenvalue of Laplacian matrix: the Laplacian matrix equals the difference between the degree matrix and the adjacency matrix. The degree matrix of a graph is a square diagonal matrix whose ii-elements equal the valency of the i-th vertex. The number of subsets of the graph that are mutually not connected equals to the number of zero eigenvalues of the Laplacian matrix. Hence for a connected graph the Laplacian matrix has exactly one zero eigenvalue. [43] Algebraic connectivity is the second smallest eigenvalue of the Laplacian matrix and always positive. It is well-known to characterize how well-connected a graph is. [45] Graph energy equals to the sum of all eigenvalues of the adjacency matrix [45].

Construction of representative volume elements
To examine the respective importance of the various descriptors introduced in the previous section, we studied a large number of representative volume elements (RVEs) with fiber networks and examined how the descriptors of these fiber networks correlated with their mechanical properties. To minimize finite size effects in our RVEs, we developed an algorithm mimicking a periodic material structure. This algorithm placed around the main RVE some additional identical ghost RVEs ( Fig. 1) that were taken into consideration only when in our main RVE descriptors had to be evaluated at nodes or fibers connected across the boundaries of our main RVE. In general, we focused on the construction of isotropic networks.
To study the correlation between network descriptors and network properties it is necessary to construct networks whose geometry complies with certain predefined values or distributions of the descriptors of interest. To construct such networks we adopted the procedure described in the following subsections.

Generation of initial network
To construct tailor-made RVEs, we started from fiber networks with a random initial configuration. To speed up the stochastic optimization procedure used subsequently to transform this initial configuration into one exhibiting the desired values and distributions of the descriptors of interest, we ensured that already in the initial configuration some descriptors of the network agreed with the desired ones, namely, its connectedness, valency distribution, maximum fiber length and node density. We chose these descriptors because it is easy to directly construct networks where these descriptors take on desired values and distributions so that it would be a waste of computational time to include these descriptors in the stochastic optimization algorithm described below. Rather this algorithms can then easily be constrained in such a way that it maintains over all iterations the initial values and distributions for these descriptors. To prescribe the above mentioned descriptors already in the otherwise random initial configuration we applied the following procedure.

Node density, connectedness, maximal fiber length
First, we defined a cubic domain of edge length L RVE for our RVE to be constructed. To keep finite size effects acceptably small, we generally imposed without loss of generality L RVE /3 as a hard constraint for the maximum fiber length (l max ). To keep this constraint already in the initial configuration and ensure a largely homogeneous initial configu- In the next step, we computed the number of nodes N node = ρ node L 3 RVE corresponding to the desired node density and distributed them as equally as possible among the subdomains to ensure a largely homogeneous initial configuration ( Fig. 2(b)). Note that an exactly equal distribution was not always possible given the predefined node density and number of subdomains, which could be achieved only if in certain cases some subdomains were assigned one node more than others.
In the third step, we connected the nodes within each subdomain by a random polygon chain of initial fibers ( Fig. 2(c)), which automatically complied with the given maximum length due to the above chosen edge length of the subdomains.
In the fourth step, we interconnected in each subdomain one node with one node in a neighboring subdomain such that a connected network was achieved ( Fig. 2(d)) and the added fibers did not violate the maximum length criterion.

Valency distribution
Let the valency of node i be v i and p v (v i ) be the desired probability distribution. The maximum valency allowed is v max , then the number of fibers in the network should be Having computed this quantity, fibers were randomly added between pairs of nodes in the network, checking each time whether the maximum length criterion was satisfied and also whether the respective addition of a fiber helped to bring the actual valency distribution closer to the desired one p v (v i ). Only if so, fibers were actually added. Otherwise another connection between a random pair of nodes was examined. To account for the periodic boundary conditions of our RVE, we allowed also connections from a node inside our RVE to a node in one of the neighboring ghost RVEs (Fig. 1). If such connections were established, the node in the ghost RVE was effectively replaced for all further considerations of the network connectivity by its periodic counterpart within the RVE.

Simulated annealing
In the previous sections, we pointed out how to construct RVEs of connected fiber networks with a maximal fiber length and a prescribed node density and valency distribution. Here we point out how to transform these RVEs into a final configuration in which also the fiber length distribution and the direction cosine distribution match desired target distributions. To this end, we use the simulated annealing (SA) method [34,4]. SA is a (global) stochastic optimization method. The idea of the approach is to define a cost function that becomes minimal if all N d descriptors of interest of the RVE match their target values or distributions. E k is the cost function (energy function) of the k-th descriptor and w k its scalar weight. For reasons discussed in more detail below, it turned out to be sufficient for our purposes to consider only two descriptors in (7), namely, the fiber length distribution and the direction cosine distribution. Following [4], we defined the cost function of the k-th descriptor via the Cramer-von Mises test, that is, Here we assume that the k-th descriptor is a probability density distribution with N k realizations x ki across the whole network, ordered such that x k1 < x k2 < x k3 < …. For example, if the k-th descriptor is the fiber length distribution, N k is the number of fibers in the network, and x ki is the length of the i-th fiber. c k is the cumulative target distribution of the k-th descriptor.
Once a cost function E has been defined that becomes minimal if the network has reached a state where the descriptors match their targets, one uses in SA a procedure of random evolution steps with the aim to decrease E down to its (ideally global) minimum. In our work we used two different types of evolution steps described in the following.
Step type 1 applies a random displacement vector (with a constraint for its maximal absolute value) to a randomly selected vertex (node) in the network (Fig. 3(b)).
Step type 2 deletes a pair of fibers and then adds a new one in such a way that the valency of all the affected nodes remains unaltered. For the configuration in Figs. 3(a), the two alternative ways to  accomplish such a step (once a pair of fibers to be deleted has been chosen) are depicted in Fig. 3(c) and (d). For step type 2 we generally checked whether it altered the connectedness of the network and admitted only such steps of type 2 that did not.
It is important to note that neither of the above two types of random steps changes the connectedness nor the valency nor the node density of the network (which we ensured to match their respective targets already in the initial configuration). However, other descriptors such as length and direction cosine between neighboring fibers change. After each random step, we computed ΔE, the change of the cost function E. Using a Metropolis algorithm as already [34], the probability to accept the step was computed in the i-th evolution step as Here T i is a temperature-like parameter. If the energy of a proposed random step decreases or at least not increases the cost function E, it is always accepted in the Metropolis algorithm. However, to avoid getting trapped in local minima, it is essential to endow the algorithm with the ability to perform at least with some likelihood also steps into an energetically unfavorable directions. This ability is ensured by the temperature parameter T i . The higher T i the more likely it is that the Metropolis algorithm every so often performs also energetically unfavorable steps. Typically, one starts with high T i to give the algorithm the chance to explore large parts of the state space and decreases T i then with an increasing number of random evolution steps [47]. Different approaches are used in the literature for this so-called cooling schedule. We applied the power function with the initial temperature T 0 and d r being the temperature decay-rate. SA is stopped if either a predefined maximal number of random evolution steps i max is reached or the cost function is lower than a predefined threshold E target . We used a combination of both conditions (for details see Fig. A.12).

Definition of RVE set used in our study
As discussed already in [4], fiber length, valency and direction cosine distributions have been reported at several occasions to affect the mechanical properties of fiber networks. It appears mechanically plausible that also the number of fibers plays an important role. For a given valency distribution this number can be expressed equivalently by the node density. The primary objective of this study is identifying the descriptors governing the mechanical behavior of fiber networks. Therefore, we decided to construct for our study a large set of RVEs which sampled specifically variations with respect to the above four descriptors, namely, valency distribution, fiber length distribution, direction cosine distribution and number of fibers in the RVE.
To this end, we discretized the space of possible valency, fiber length, and direction cosine distributions by a set of trial distributions illustrated in Fig. 4. The algebraic formulae for these distributions denoted by V i , L j and C k with i = 1, 2, …, 5, j = 1, 2, …, 6 and k = 1, 2, …, 5 are specified in Appendix B. The set of trial functions was defined heuristically, though in a manner that allowed us to sample with a relatively small number of trial functions a large part of the possible shapes these functions typically take on in real physical systems. To this end, we included for each descriptor at least a constant, two linear and a quadratic trial function. Moreover, we allowed an independent variation of the number of fibers F q with q = 1, 2 between the two specific values F 1 = 4000 and F 2 = 8000. Note that we numbered the considered distributions such that the mean value of the respective descriptor increases monotonically with their subscript. For instance, L 1 refers to the distribution with the smallest average length and L 6 to the distribution with largest one.
Our definition of trial functions allowed us to define 300 types of different RVE, each characterized by a specific choice of the descriptors V i ,L j ,C k ,F q . These different types of RVE are referred to in the following by the abbreviations V 1 L 1 C 1 F 1 ,V 2 L 1 C 1 F 1 ,…, V 5 L 6 C 5 F 2 , respectively. It is worth mentioning that SA did not converge satisfactorily for the combination of V 3 , V 4 and V 5 with C 1 and C 4 , suggesting that these types of

Table 1
Possible combinations of V i , L j and C k . RVE cannot exist due to a violation of intricate, problem-intrinsic geometric constraints. Similarly, RVE with fiber length distributions of L 1 and L 6 could be constructed only for V 1 in combination with all considered direction cosine distributions except C 1 . RVE types that could be constructed with F 1 could always be constructed also with F 2 . Hence, among the theoretically possible 300 different RVE types only 168 could be constructed, which are summarized in Table 1. For each RVE type, 15 different random realizations were generated by the procedure described in Section 3. Thereby, we used the SA parameters from Table 2 for the construction of all the 2520 RVEs considered in our study. Fig. 5 illustrates the construction of one specific such RVE (of type V 5 L 2 C 2 F 1 ). Starting from a random initial configuration, SA yields an RVE where the descriptors of interest closely match the prescribed target distributions.

Descriptors governing mechanical properties
To understand the influence of the descriptors of the network structure on the network's mechanical properties, we discretized the network RVE generated as described in Section 3.3 by finite beam elements, to which we assigned a circular cross section with a diameter of L RVE /50, a Young's modulus of E f = 79 GPa and a Poisson's ratio of ν f = 0.44. These values are deliberately chosen to resemble those of nanoporous gold as a possible application case. Yet, the exact choice of these parameters is largely irrelevant for the following discussion as long as they enable in principle all relevant deformation modes of the network in a physically realistic range, which we confirmed to be the case. Fibers were discretized using beam finite elements based on the Timoshenko beam theory. At the intersection points of fibers, both translations and rotations were coupled, i.e., rigid joints were assumed. Applying periodic boundary conditions to our RVE according to [48], we computed their (effective, homogenized) Young's modulus E and Poisson's ratio ν by the method suggested in [28].
To investigate the effect of valency distribution, fiber length distribution, direction cosine distribution and the number of fibers on mechanical properties, we divided the 2520 RVEs generated according to Section 3.3 into categories. V i L j C k F q denotes a category of an RVE where all the four descriptors match their respective target distributions. That is, the RVE category includes the 300 RVE types V 1 L 1 C 1 F 1 ,V 2 L 1 C 1 F 1 ,…, V 5 L 6 C 5 F 2 . Analogously, by leaving out one or several of the descriptors in the list, we denote an RVE category where only the listed descriptors were fixed to certain target distributions and the others were allowed to vary freely among all distribution functions included in our study. For example, L j C k denotes a category of RVEs where the length distribution matches L j and the direction cosine distribution C k , whereas all the other descriptors may be arbitrary (within the set of distributions and values considered herein). As there are altogether 30 different RVE types of the category L j C k , denoted by L 1 C 1 , L 1 C 2 , L 1 C 3 , …, L 6 C 5 , our study includes for each of these RVE types 2520/30 = 84 realizations (in case all combinations of remaining descriptors be available for the subsets of this category). In the extreme case that all descriptors are allowed to vary freely among the distributions included in this study, we denote the respective RVE category (or type) by '-', including all 2520 RVE generated in our study. As illustrated also in Fig. 6, there are altogether  3 MPa of the RVE type V 2 L 3 C 2 F 1 ); (b) fixing more and more descriptors, the relative variance (σ 2 ) of the Young's modulus in the resulting RVE categories (boxes) decreases by the percentage indicated next to the lines that connect an RVE category with one that results from fixing a further descriptor. 16 RVE categories, which are listed on the very left of Fig. 6. For each RVE category, we computed for all RVE types belonging to it the mean values of Young's modulus and Poisson's ratio as well as the variances of these quantities (normalized by their respective mean value) (σ 2 ). For each RVE category the average variance across all the RVE types belonging to it is plotted in Fig. 6. Moreover, we computed in each RVE category for each RVE type belonging to it for all samples the maximal absolute values of the relative deviations of Young's modulus and  Poisson's ratio from their mean values. This quantity Δ max was averaged for each RVE category and plotted in Fig. 6.
In general, the larger σ 2 and Δ max for a category, the less do the descriptors fixed for this specific RVE category determine the mechanical properties of a network. By contrast, in the theoretical extreme case that σ 2 and Δ max are zero, the descriptors fixed in the related RVE category fully determine the mechanical properties of the fiber network (at least in the theoretical case that for each RVE type an infinite number of realizations is included in the study). Fig. 6 reveals that for the RVEs in the category '-', the relative variance of Young's modulus across all the subset is as big as 188.9%. By contrast, if all the four descriptors valency distribution, fiber length distribution, direction cosine distribution and the number of fibers are fixed, that is, in the category V i L j C k F q , the relative variance of the Young's modulus is as small as 0.5%. For Poisson's ratio the relative variance in the category V i L j C k F q takes on the similarly low value of 0.3%. However, for Poisson's ratio the relative variance even in the category of RVE without any geometric constraints is only 2.6%. Together, this leads to the following main conclusions: Poisson's ratio of (initially stress-free) isotropic random fiber network generally exhibits a relatively low statistical variance with a mean value of around 0.25. By contrast Young's modulus strongly depends on the geometry of the fiber network. Nevertheless, the four descriptors valency distribution, fiber length distribution, direction cosine distribution and the number of fibers together suppress around 99.5% the variance of the linear-elastic mechanical properties of fiber networks. Fig. 7(a) analyzes specifically the effect of F q by plotting for a representative selection of RVE types on the Young's modulus normalized by the number of fibers. Apparently, this ratio is (except for minor supposedly mainly statistical deviations) constant for each RVE type. In other words, if the other descriptors are fixed, Young's modulus is directly proportional to the number of fibers. This simple law can also be understood analytically. Imagine an RVE with a specific size and F q fibers. The network in this RVE can be modeled as an elastic spring S. Increasing the number of fibers, for example, by a factor of two can be constructed as process where we add another random fiber network of the same type into the RVE so that the RVE volume is then occupied by two interpenetrating random fiber networks. As these are of the same type, they exhibit the same elastic properties. That is, the RVE harbours then two elastic springs of the kind of S that act in parallel, increasing Young's modulus by a factor of two. By contrast, as both elastic systems of type S exhibit the same deformation behavior, Poisson's ratio is not expected to change if the number of fibers increase by a factor of two. This explains the very low sensitivity of Poisson's ratio to variations of number of fibers observed in our study. Of course, this illustrative explanation is applicable to any scaling factor for F q , yielding directly the simple linear relation between F q and Young's modulus observed in Fig. 7(a).
For the RVE category F q , that is, the one with fixed number of fibers, Fig. 7(b) illustrates how fixing more and more descriptors suppresses in a step-wise manner nearly the whole variance of Young's modulus. Apparently, the mechanically most important single descriptor in this process is the valency distribution, which alone reduces the variance of the mechanical properties by nearly 90% when fixed. In fact, this interpretation qualitatively still holds if we examine in Fig. 7(b) arbitrary transitions from one RVE category to another by fixing one additional descriptor. In every case, fixing the valency distribution reduces Fig. 10. Comparison of (a) normalized variance σ 2 and (b) the maximal absolute value of the relative deviation Δ max for the first three geometric moment invariants I 1 , I 2 and I 3 of the RVEs in the two categories '-' and V i L j C k F q . Fig. 11. Comparison of (a) spatially average normalized variance σ 2 and (b) the spatially averaged maximal absolute value of the relative deviation Δ max of c pore (r) , c r− node (r) and c r− fiber (r) for RVEs of the categories '-' and V i L j C k F q . the variance more than fixing any other descriptor. This justifies the conclusion that the valency distribution is generally more important than fiber length distribution and direction cosine distribution.
In order to better understand the influence of the remaining descriptors on the mechanical properties of fiber networks, Young's modulus is plotted in Fig. 8(a)-(c) for a number of representative RVE types. Fig. 8 illustrates that changes in the valency distribution are associated with changes of Young's modulus by more than one order of magnitude whereas changes in the fiber length and direction cosine distributions typically have a much smaller effect. As the indices of the descriptor distributions increase with their mean values, Fig. 8(a)-(c) demonstrate for a representative selection of RVEs that Young's modulus increases with the mean valency V mean but decreases with mean fibers length L mean and mean direction cosine C mean . To quantify this dependence, Fig. 8(d)-(f) uses a special post-processing of the results of our computational study. Thereby, we first specified a certain descriptor of interest (either valency, fiber length or direction cosine distribution) for which we studied its relation to Young's modulus. To do so, we selected from the 16 RVE categories listed in Fig. 6 the ones where the respective descriptor of interest was fixed. For example, if valency distribution is selected as the descriptor of interest, these are the eight For each of these categories, we defined subcategories containing the RVE types where the descriptor of interest was varies across all its possible distributions while all the other descriptors were kept constant at specific values. For example, the category V i F q was divided into a subcategory containing the RVE types varying the descriptor of interest V i and keeping F q constant at F 1 ), and another subcategory with the RVE types varying the descriptor of interest V i and keeping F q constant at F 2 ). For each subcategory we computed for all RVE realizations included in our study Young's modulus normalized by the mean Young's modulus within that category E ref as well as the mean value of the descriptor of interest. We excluded subcategories where not all RVE types were available because they could not be constructed (see Table 1). The results are depicted in Fig. 8(d)-(f), where the descriptor of interest was chosen to be the valency distribution, fiber length distribution and direction cosine distribution, respectively. Interestingly, Fig. 8(d) reveals a quadratic relation between the Young's modulus and the mean value of the valency distribution. This relation is remarkably linear between the Young's modulus and the mean values of the fiber length distribution and the direction cosine distribution (Fig. 8(e) and (f), respectively). The box plots illustrate the statistical deviations from the mentioned relationship in each of the above defined subcategories. The apparently very small size of the boxes reveals that the defined function almost completely characterizes the relation between the considered descriptors and Young's modulus. This observation is remarkable because it means that in discussing the relation between microstructure and mechanical properties of the RVEs, one can focus not only on four key descriptors in the form of probability distributions but in fact largely on mean values. In other words, the linear mechanical properties of the RVEs are largely determined by the number of fibers (per volume) and the scalar mean values of the valency, fiber length and directioncosine distribution. From Fig. 8(d)-(f) it is also apparent that, as discussed already above, mean valency is of much higher importance than mean fiber length and mean direction cosine with respect to the resulting material stiffness. Fig. 8(e) and (f) cover a large part of the range within which fiber length and direction cosine can be varied in a physically and geometrically reasonable way. Yet, they display variations of Young's modulus only by approximately a factor two and three, respectively. By contrast, in Fig. 8  (d) we observe a factor of almost 70. Additionally, it is worth mentioning that the valency range up to slightly above five is well-suited for random fiber networks. However, for ordered ligament systems with hexagonal cells or tetrahedral cells even higher mean valencies may appear (for example six for a uniform cubic mesh and twelve for a uniform tetrahedral mesh), which underlines even more impact on Young's modulus that can in principle be realized in networks of ligaments by varying the mean valency.
As revealed already by Fig. 6, Poisson's ratio is nearly constant across the different RVE types, which is why we omit a detailed discussion of its minor dependencies on the different descriptors. We note, however, that we have shown already above that Poisson's ratio can be expected to be (nearly) independent on the number of fibers F q . Given that F q is generally responsible for a large part of the variations of the mechanical properties observed in this study, the insensitivity of Poisson's ratio to it may explain why Poisson's ratio generally exhibits only relatively small variations in the RVE studied herein.
According to [49][50][51], the ratio of the fiber radius to the fiber length plays an important role for the mechanical properties and deformation of random fiber networks. In networks with higher values of this ratio, stretch is the dominant deformation mode which causes the network to deform affinely. In contrast, if the ratio is small, bending is the predominant deformation mode and the network deforms non-affinely. In order to analyse the dependence of our conclusions in this section on geometrical properties of the fibers such as their slenderness ratio, we repeated the studies discussed in this section with networks with a fiber diameter of L RVE /500. The results are presented in Appendix C. They exhibit some minor quantitative differences to the ones obtained with a fiber diameter of L RVE /50 but are qualitatively similar, underlining thus the robustness of our conclusions.

Relation between morphological descriptors and graph descriptors
In Section 4.1 we demonstrated that a set of four geometric descriptors nearly fully determines the mechanical properties of random fiber networks. However, such networks can also be interpreted as mathematical graphs, the nodes forming the vertices of the graph and the fibers its edges. Mathematical graph theory provides a host of descriptors for graphs. It is instructive to study the relation between such graph descriptors and the four geometric descriptors on which Section 4.1 focuses. In this section we examine the graph descriptors introduced in Section 2. As classical graph theory works with finite rather than infinite graphs, most descriptors from graph theory cannot deal with the assumed periodicity of the considered RVEs. Therefore, we evaluated them only within the RVE, ignoring the fibers cutting through the boundaries of the RVE.
In order to investigate the dispersion of graph descriptors, we evaluated the aforementioned statistical quantities σ 2 and Δ max also for a host of graph descriptors and compared the results for the two extreme cases, the RVE category '-' and the RVE category V i L j C k F q . The results are illustrated in Fig. 9. The ones related to V i L j C k F q are the mean values of the 168 types RVE belonging to the category. Fig. 9 reveals that fixing the four descriptors valency distribution, fiber length distribution, direction cosine distribution and number of fibers reduces σ 2 and Δ max for all graph descriptors substantially. Except for the clustering coefficient, all graph descriptors seem to be almost fully determined in an implicit manner by fixing the four dominant morphological descriptors, underlining once more their key role in characterizing properties of the fiber network.

Role of other morphological descriptors
We have established now that fiber length distribution, valency distribution, direction cosine distribution and number of fibers are the four morphological descriptors that nearly fully determine the mechanical properties (and indeed also graph properties) of fiber networks. In this section we analyze, to which extent they determine also the other morphological descriptors introduced in Section 2. The first three geometric moment invariants, I 1 , I 2 and I 3 [42] are computed for the RVE types within the two RVE categories'-' and V i L j C k F q . Analogously to the previous sections, the mean values of the normalized variance σ 2 and the maximal absolute value of the relative deviation Δ max averaged over all RVE types in these two categories are plotted in Fig. 10. Apparently, fixing the four above identified morphological main descriptors determines also the three first geometrical moment invariants nearly completely. Higher order moments exhibit very small values compared to the first three ones. They are typically considered as of minor importance compared to the first three ones and thus skipped here.
Furthermore, the pore-size distribution function p pore (r), the radial distribution function of the nodes p r− node (r) and the radial distribution function of the fibers p r− fiber (r) have been studied in this work. To determine the pore-size distribution function, we chose 1000 random points inside the void-phase of the RVE and computed their respective distance to the closest fiber. Then the probability distribution of distances was approximated in a discrete manner by dividing the occurring interval of 0 < r/L RVE < 0.12 into 24 bins with equal width. A similar calculation was performed for p r− node (r) and p r− fiber (r) by dividing the occurring interval 0 < r/L RVE < 2 into 25 equal-sized bins. For the computation of p r− fiber (r), fibers were divided into segments with a maximum length of 5%L RVE . In order to compute statistical quantities comparable to σ 2 and Δ max above, we first converted the probability density functions (PDFs) p pore (r), p r− node (r) and p r− fiber (r) into associated cumulative (CDFs) distribution functions c pore (r), c r− node (r) and c r− fiber (r).
For these, we calculated σ 2 and Δ max separately in each bin used for the discretization of the distributions and introduced in general for the characterization of a CDF c(r) the averaged quantities Here, L c = 0.12L RVE for c(r) = c pore (r) and L c = 2L RVE for c(r) = c r− node (r) and c(r) = c r− fiber (r). Fig. 11 compares σ 2 and Δ max for the mentioned distribution functions for the RVEs in the categories '-' and V i L j C k F q . Apparently also p pore (r), p r− node (r) and p r− fiber (r) are at least to a large extent implicitly determined by fixing the four morphological main descriptors identified above.

Conclusions
Studies of the mechanical properties of heterogeneous random media have attracted substantial interest over the last two decades. A key question in this area is which descriptors of the microstructure of such media determine to which extend their mechanical properties. Specifically for networks of fibers or ligaments, this has remained unclear so far. Such networks play an important role in biophysics and soft matter physics but also in materials research, in particular with respect to nanoporous metals [52][53][54][55][56]. The microstructure of the latter differs from the fiber networks usually studied in biophysics by the fact that it is formed by a network of ligaments with a considerable thickness and variation in shape. Yet both fiber networks in biophysics and ligament networks as occurring in nanoporous metals appear to share important common properties. For example, in biophysics, it has been suggested that fiber length distribution, valency distribution and direction cosine distribution play important roles in the mechanical behavior of collagen fiber networks [4]. And also for nanoporous metals it has been found that characteristics such as node valency or ligament length have important effects on the macroscopic mechanical properties [57][58][59][60]. Yet, so far it remains unclear to which extend these findings in separate areas of materials science are related to each other and which other descriptors are important determinants of the mechanical properties.
To answer these questions, we performed a large computational study including more than 2500 RVEs. Using simulated annealing [33], we constructed these RVEs sampling a large part of the physically reasonable configuration space. As expected, we observed in finite element models a large variance of Young's modulus and Poisson's ratio among these 2520 RVEs. However, we also observed that fixing four morphological key descriptors to specific values or probability distributions reduced this variance by more than 99%. These four key descriptors are the number of fibers per volume (i.e., the fiber density), the node valency distribution, the fiber length distribution, and the direction cosine distribution. We demonstrated that the four key descriptors did not only largely determine the mechanical properties but also a host of other morphological or graph-theory-based descriptors. We thus conclude that the number of fibers per volume, the node valency distribution, the fiber length distribution, and the direction cosine distribution are the key properties of networks of (thin) fibers both with respect to mechanics and geometry.
Remarkably, we observed that the number of fibers per volume is (nearly) linearly related to Young's modulus of fiber networks, which we could also explain by a simple theoretical analysis. Interestingly, we observed that also the effect of valency, fiber length and direction cosine distribution on Young's modulus can be captured to a very large extent by a simple linear or quadratic relationship between the mean values of these distributions and Young's modulus. Thereby stiffness increases with mean valency and decreases with mean fiber length and mean direction cosine. This finding is in excellent agreement with [61], where based on a different approach already previously a remarkably simple ascending relationship between mean valency and stiffness was reported specifically for nanoporous metals.
Given that for valency, fiber length and direction cosine distribution apparently mainly the mean values matter for the stiffness of the network, one can summarized our study as follows: the mechanical properties of fiber networks are largely determined by four scalar descriptors, which are the number of fibers per volume and the mean valency, mean fiber length and mean direction cosine. Among these four scalar descriptors the former two are by far dominant whereas the latter two are modulators of minor importance. Generally, we found that the dependence of Young's modulus on the four key descriptors was much stronger than the dependence of Poisson's ratio.
The identification of four simple scalar descriptors for characterizing the mechanical properties of fiber networks can be expected to be an important step to understand the relation between microstructure and macroscopic mechanical properties of materials consisting of random fiber networks. It should, however, be kept in mind that this study has also certain limitations. First, it focuses on linear mechanical properties (Young's modulus and Poisson's ratio) of isotropic fiber networks only. A generalization to nonlinear anisotropic materials is a natural next step. Since the macroscopic material behavior is strongly influenced by microstructural defects [62][63][64], large deformation, plasticity and failure of random fiber networks is a promising area of future research. Second, our study largely relies on a computational approach. A promising next step is underpinning its main findings by a careful theoretical analysis unraveling the rationale behind the role of the different descriptors.

Data availability
The raw/processed data required to reproduce these findings cannot be shared at this time as the data also forms part of an ongoing study.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
The probability distribution functions (PDFs) for the valency, fiber length and direction cosine used in this article are defined by algebraic equations in the following.

B.1. Valency distributions
where v (v = 1, 2, …, 6) represents the node valency. where c ( − 1⩽c ≤ 1) is the cosine of the angle between a pair of fibers adjacent to the same network node.  4e3 MPa of the RVE type V 3 L 3 C 3 F 1 ); (b) fixing more and more descriptors, the relative variance (σ 2 ) of Young's modulus in the resulting RVE categories (boxes) decreases by a percentage indicated next to the lines connecting an RVE category with one resulting from fixing one more descriptor, respectively.