Strain Measures and Energies for Crimped Fibres and Novel Analytical Expressions for Fibre Populations: Ingredients for Structural Fibre Network Models

This paper deals with the structural modelling of fibre networks with a focus on the description of populations of initially crimped fibres. It presents a systematic approach of introducing appropriate strain measures for single fibres based on a deformation decomposition and by transferring knowledge from the field of elastoplasticity. On this basis, for example, the often used Biot-type fibre strain measures λ−λw\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\lambda-\lambda_{\mathrm{w}}$\end{document} and λ/λw−1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\lambda /\lambda _{\mathrm{w}}-1$\end{document}, with stretch λ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\lambda $\end{document} and waviness λw\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\lambda _{\mathrm{w}}$\end{document}, are consistently assigned to different classes of material strain (“additive”) or intermediate strain (“multiplicative”) descriptions, respectively. We review in this work different fibre strain energies based on the different stain measures and present extensive comparisons on the physical implications and the results on the fibre population and network scale. These investigations also include formulations with a Hencky-type energy based on a logarithmic strain. Furthermore, we present novel analytical expressions for fibre populations that make the evaluation of integral expression superfluous and thus lead to a significant reduction in computational time.

elastin and collagen fibres. The so-called full network models are particularly well suited for the multi-scale description and modelling of such materials within the framework of continuum mechanics. Such models have been formulated in the context of rubber-like materials by, for example, Treloar [121], Treloar and Riding [122], Wu and van der Giessen [130] and Miehe et al. [88]. For biological materials, such network models were proposed in the fundamental works of Lanir [72,73] and have been employed for the modelling of, for instance, tendons and ligaments by Hurschler et al. [63], articular cartilage by Federico and Herzog [39] and Ateshian et al. [5], blood vessels by Alastrué et al. [2] and skeletal muscles by Bleiler et al. [16]. Moreover, the conceptually similar microplane models by Bazant and Oh [10] and Carol et al. [24] have been applied by Caner and Carol [22] to the modelling of blood vessels and by Caner et al. [23] to the intervertebral disc. In brief, network models proceed from a parametrisation of space directions by means of spherical coordinates and the calculation of the effective network response as a result of the mechanical properties of the individual space directions. Network models for biological tissues, which are the focus of this paper, usually consider a space direction as a population of fibres (or fibre bundles). These populations are characterised by fibres that have different degrees of crimp in the unloaded reference configuration. Since the thin fibres mainly store elastic energy when they are stretched in their uncurled state, models usually account for the waviness of the fibres as structural information. In turn, the structural properties of fibre populations are described in a statistical manner via so-called waviness distribution functions, as done, for instance, by Lake and Armeniades [71], Decraemer et al. [29], Lanir [73], Horowitz et al. [61], Hurschler et al. [63] or Sacks [109]. Further, in order to quantify the stretch of single fibres beyond the waviness, appropriate strain measures have to be employed. Different strain formulations have been used in previous studies such as, for example, λ − λ w [e.g., 21,27,71] and λ/λ w − 1 [e.g., 29,63,73], with stretch λ and waviness λ w , However, despite the numerous applications of network models, the different meaning of those two, their influence on the results and further aspects related to fibre populations have not been addressed in detail so far.
Thus, we focus in this contribution on some open theoretical topics and structure the paper as follows: The basics on network models are briefly repeated in Sect. 2. In Sect. 3, we describe the modelling of populations (bundles) of unidirectionally-oriented fibres with different initial degree of crimp and address the requirements for appropriate waviness distributions functions. Subsequently, after looking at fibre networks and fibre populations, we consequently move down the length scales and focus in Sect. 4 on the structural modelling of individual crimped fibres. Therein, we propose, to the best of our knowledge, for the first time a systematic way of introducing strain measures for crimped fibres. To do so, we use a deformation decomposition and knowledge from the field of elastoplasticity. The procedure results in the formulation of two strain measure families, one referred to as additive and one as multiplicative. Subsequently, we discuss commonly employed fibre energy functions and show the different inherent physical meanings. In Sect. 5, we present novel analytical expression for fibre populations that make the numerically expensive solution of integral expressions superfluous. Finally, examples and related discussions are provided in Sect. 6. In particular, we compare the different fibre energy formulations on the level of fibre populations and the entire fibre network. Section 7 draws some concluding remarks.

Network Models and Effective Quantities
First, we briefly introduce the basics of network models. This provides the context for the following investigations on populations of unidirectionally oriented fibres.

Basics of Network Models
Starting point for the network models is the idealisation of the fibrous (or chain-like) microstructure of a material by assuming that the fibres connect the origin of a unit sphere with some point on the associated sphere surface, such that the orientations of the fibres coincide with the ones in the real microstructure, see Fig. 1. The fact that the unit sphere describes the microstructure of a material motivates to call it microsphere, see, for example, Miehe et al. [88] or Alastrué et al. [2]. Subsequently, it proves useful to introduce a unit vector r 0 ∈ R 3 that connects the centre of the microsphere with a point on its surface, denoted by , and that is parametrised using spherical coordinates, leading to Therein, ∈ [0, π) denotes the polar angle, ∈ [0, 2π) is the azimuthal angle, and the e s i -coordinate system (i = 1, 2, 3) represents a laboratory frame of reference, see Fig. 2a. From Eq. (1) follows that the Euclidian norm of r 0 is |r 0 | = √ r 0 · r 0 = 1. The vector r 0 is regarded as a cylinder-like segment of the microsphere that is occupied by a population of unidirectionally oriented fibres (also called fibre bundle), see Fig. 2b. Subsequently, it is essential to describe the micro-kinematics of the fibre population (hence, the vector r 0 ). This is done by means of a general mapping operator F r : r 0 → r = F r (r 0 ) that provides the vector r in the actual configuration. In particular, it is convenient to postulate the existence of a linear map F r ( , ) : r 0 → r = F r r 0 , where the second-order tensor F r ( , ) describes the micro-deformation of the referential vector r 0 , see, for example, Fried [44] or Ehret [32] After having defined the basic kinematics of the network model, the mechanical properties of a fibre population are described. Since we focus on the purely hyperelastic case, the mechanical behaviour is fully described by assigning a strain-energy function W r to each vector r 0 ( , ). The potential-like strain energy is a frame-indifferent function of the micro-deformation F r , such that W r (QF r ) = W r (F r ) for all proper orthogonal tensors Q ∈ SO (3). A convenient (but not necessary) way for ensuring frame-indifference is to formulate W r in terms of invariant deformation measures such as where the micro-stretch λ r describes the longitudinal stretch, the area stretch υ r quantifies the cross-sectional deformation and the micro-Jacobian J r the volumetric deformation of the cylinder-like segment, see, for example, Miehe et al. [88], Carol et al. [24] or Govindjee et al. [49]. Network models for collageneous tissues usually assume purely stretch-based energies W r (λ r ). However, υ r -based energies are common in the field of rubber elasticity and lead, for example, to so-called tube-like topological constraints as formulated by Heinrich and Straube [55] and Edwards and Vilgis [31]. Moreover, J r -dependent volumetric energy contributions are classical ingredients in rubber elasticity, see, for example, Wall and Flory [125]. They can be interpreted as inter-chain repulsive forces, see Bischoff et al. [15]. As we focus in this work on collageneous networks, we proceed with stretch-dependent strain energies W r (λ r ).
The key idea of the network model is that the orientations of the fibres in the idealised microsphere shall be identical to the ones in the real microstructure. Hence, there is a certain probability of finding the end point of a fibre on a specific point of the microsphere surface Ω. By further postulating that the number of fibres is large enough to assume a continuous distribution of fibres in the microsphere, the orientation of the fibres can be described in a statistical sense by means of an orientation density function (ODF) p (r 0 ) = p ( , ). The ODF describes the one-point probability that the vector r 0 points to an infinitesimal area element dΩ = sin[ ] d d ∈ Ω of the microsphere surface Ω. Considering the spherical coordinates and the periodicity of the angles and , the ODF has to satisfy the conditions p (r 0 ) = p (−r 0 ) and p ( , ) = p ( + nπ, ) = p ( , + 2nπ) with n ∈ N . These conditions suggest that the field of directional (also called spherical) statistics has to be consulted and the ODF has to be defined in terms of bivariate spherical probability functions, see, for example, Mardia and Jupp [82]. For completeness, it is further remarked that p can also be chosen as a Dirac Delta function with a certain number of masses, which means that the continuous fibre distribution degenerates to a discrete number of fibres. This procedure would lead to classical discrete chain models that are used in polymer mechanics, such as the three-chain models by James and Guth [65] and Wang and Guth [126], the four-chain model by Flory and Rehner [40] and Treloar [120], or the eight-chain models by Arruda and Boyce [4] and Kroon [69], see also Beatty [11].

Effective Affine Network Response
Having defined the micro-energy W r (F r ), the averaged stored energy in the microsphere is obtained as Ω = W r (F r ) Ω , where (·) is a continuous averaging operator over the microsphere surface defined as It is consequent to formulate a normalisation condition for p such that 1 Ω = Ω p ( , ) dΩ = 1. In classical homogenisation treatments the micro-deformations F r are such that Ω is minimised and the effective energyW is defined in terms of a minimum energy principle, see, for example, Hill [57] or Ponte Castañeda and Suquet [101]. We do not apply such approaches in this work, but, for completeness, we provide some further details on variational energy minimisation methods that lead to nonaffine micro-deformation in Appendix A. Here, we instead proceed from an affinity assumption, which means that each vector r 0 experiences the macroscopic deformation. This leads to the affine map With this, the micro-deformation measures from Eqs. (3) can be reformulated tō Note that r,λ r andν r are still dependent on and . With prescription (5), the effective network energy is directly given asW by means of a simple averaging of the fibre energies over the microsphere surface Ω. Therein,W V denotes the Voigt estimate, which arises from the assumption of uniform microdeformations. It is well known that the Voigt estimate is a rigorous upper bound for the effective behaviour, see Ogden [96], but serves as a reasonable prediction and easy-to-use tool if the affinity assumption is, for example, justified by experimental observations. In the remainder of this work, we focus on strain energy formulations W r for populations of unidirectionally-oriented fibres. Therefore, the following investigations are independent of the specific choice between a minimisation principle (described in Appendix A) and the affinity assumption. Thus, it is sufficient to proceed with a general mapping function F : r 0 → r = F r 0 and a stretch λ = |r| = |F r 0 |, such that W r = W r (λ). The results presented in the following can be transferred directly, for example, to the affine model by replacing F withF from Eq. (5) and λ withλ r from Eq. (6) 1 .

Modelling of Fibre Populations via Waviness Distribution Functions
The network models introduced in the previous section rely on the assignment of a strainenergy function W r for each space orientation r 0 . In general, the theoretical framework of network models is not restricted to a specific behaviour of space orientations and any form of strain-energy functions W r can be applied. However, as illustrated in Fig. 2b, a space orientation r 0 represents a cylinder-like element with unidirectionally-oriented single fibres. In soft biological tissues one usually observes a typical J-like stress curve with gradual stiffness increase upon tensile stretch λ > 1, while there is hardly any resistance under compressive stretches λ < 1. This behaviour is attributed to the fact that the fibres-in case of soft biological tissues, these are usually collagen fibres-are not straightened in the reference configuration, but are wavy and exhibit a certain crimp, see, for instance, Wolinsky and Glagov [129], Lake and Armeniades [71] or Roy et al. [107]. This means that the initial contour length, say l f , of a fibre is longer than the fibre's initial end-to-end length, say l 0 . In turn, the different degrees of crimp of single fibres with the same end-to-end orientation r 0 entail that the fibres become straight (recruited) at different stretches λ and that the overall behaviour of the space direction r 0 becomes successively stiffer with increasing stretch. This is depicted in the diagram on the left hand side of Fig. 3.
The mechanical description of the gradual stiffening and the characteristic J-like stressstretch curve is usually done using either (i) a phenomenological or (ii) a structural approach for the strain energy W r . Thereby, a phenomenological approach proceeds from a generic strain-energy function and the subsequent calibration to experimental data. Frequently used mathematical forms are, for example, the exponential strain-energy function by Holzapfel et al. [60], W HGO = k 1 exp[k 2 (λ 2 − 1) 2 − 1]/(2k 2 ) with material parameters k 1 and k 2 , or the power function by Balzani et al. [8], W BNSH = α 1 (λ 2 − 1) α 2 with material parameters α 1 and α 2 . Many more such stretch-based or so-called J 4 -based strain-energy functions, where the fourth invariant J 4 = λ 2 is the squared stretch, can be found in the review article of Chagnon et al. [25]. Yet, this work is concerned with structural modelling approaches for the description of W r . In contrast to phenomenological approaches, structural models describe not only the observed form of the mechanical response in a mathematical way, but rather the responsible mechanisms and structures. This is done in the following subsection.

Structural Modelling of Fibre Populations
As outlined before, the reason for the stiffening of fibre populations upon tensile stretch is the sequential straightening of initially crimped fibres with the tensile resistance of straightened fibres being several magnitudes higher than in the crimped state. Thus, the deformation-in the present case the one-dimensional stretch-at which an individual fibre becomes straight is a decisive kinematical measure. It is defined as the ratio of contour to end-to-end length via commonly referred to as the waviness of a fibre. Hurschler et al. [63] called it the "straightening stretch-ratio". The same kinematical information is expressed by the reciprocal of λ w , which is usually introduced as the straightness. It makes sense that the mechanical description of a single fibre makes use of the waviness in order to distinguish between the different characteristic deformation states, as depicted in Fig. 4. Hence, in a hyperelastic setup and Fig. 4 The different states of an initially crimped fibre upon compressive and tensile stretches along its end-to-end direction. The fibre mainly undergoes bending and geometrical changes of its contour line for stretches λ ≤ λ w , while the actual lengthening of the fibre along its contour line is induced for λ > λ w when the stretch exceeds the waviness λ w for stretch-based formulations, the strain energy of a single fibre is introduced as a function of the stretch λ and should include the waviness λ w as parameter, such that Then, it is consequent to accompany the energy formulation with stochastical descriptions for the distribution of the waviness λ w in terms of a univariate waviness distribution function p w (λ w ), following the works of Lake and Armeniades [71], Decraemer et al. [29] and Lanir [73]. Mathematically the distribution p w is a probability density function and is therefore referred to in the following as PDF. The existence of a density function postulates a continuous distribution of fibres and thus represents a limit case of the idea of a discrete number of gradually recruited springs as it was proposed in the early works of Frisén et al. [45] and Viidik [124]. Alternative terms that can be found in the literature for p w are undulation or crimping distribution function. Further, Hill et al. [56] referred to it as recruitment probability distribution, similarly Fan and Sacks [36] refer to it as distribution of the gradual fibre recruitment, Avazmohammadi et al. [6] as recruitment function and Soong and Huang [116] as collagen arrival density. In addition, exceeding the threshold stretch λ w is referred to in some studies as activation of the fibres. However, some attention is needed with this term to avoid confusion in the context of biological tissues with active (e.g., contractile) behaviour, such as muscles.
The structural modelling approach for fibre populations as explained so far is thus usually based on the two basic assumptions: (i) Uniform material properties: The fibre population consists of unidirectionally-oriented, identical fibres with uniform material properties, for instance, equal tensile stiffness. (ii) Non-uniform structural properties: Individual fibres differ in their geometrical arrangement, for instance, distinct waviness values, except for their collinear end-to-end orientation and initial end-to-end length.
Hence, the effective behaviour is a result of a constitutive description for individual fibres with uniform material properties and the stochastical description of their non-uniform structures by means of a probability density function, as illustrated in Fig. 3. Finally, the resulting strain-energy function W r for a fibre population is given by the integral formulation where λ w is the lower bound of the distribution p w . More details on the bounds as well as on the necessary normalisation property of the distribution follow in the next section. It should be noted that each fibre experiences the stretch of the space orientation r 0 . This may seem like a Voigt-based affinity assumption, as applied in Eq. (9) for the effective network energy, but it instead is a direct consequence of the kinematic description of the fibres based on their end-to-end points. The first Piola-Kirchhoff stress tensor for a fibre population is obtained via a derivation of the strain-energy function W r with respect to the deformation gradient F as by making use of the chain rule and the derivative ∂ F λ = λ −1 F r 0 ⊗ r 0 = λ −1 r ⊗ r 0 . The scalar stress P r is the derivative of W r with respect to λ. Since λ occurs in the integrand and in the upper limit, P r has to be calculated by means of the Leibniz integral rule (or Leibniz's theorem) and reads With (10), it is straightforward to show that P r is associated with a nominal traction vector P r r 0 = λ −1 P r (λ)r, whereas P r r ⊥ 0 = 0, with r 0 · r ⊥ 0 = 0. This means that a stretchdependent fibre energy W r (λ) does not induce any stress contributions transverse to r 0 .
The fourth-order nominal elasticity tensor of the fibre population is obtained by means of the second derivative as where L r (λ) := ∂ λ P r (λ) = ∂ 2 λ 2 W r (λ) denotes a scalar nominal stiffness modulus that is given by Further, (·) T 23 in Eq. (12) denotes a transposition of (·) by interchanging the second and third basis vectors. So far, the fibre population, hence the energy W r , consists of mechanically identical fibres that may only vary in their initial structural arrangement. It makes sense to postulate that mechanically distinct fibres also have different structural arrangements. This means that distinct probability functions p α w and p α have to be formulated in that case, where α indicates different fibre types. Hence, this results in separate energy functions W α r for each fibre type and the overall network energyW (orW V ) is the result of all fibre populations. In this work, we will focus on mechanically identical fibres. Further, Eqs. (11) and (13) account for general fibre energies W f and, thus, also to such that are not zero in the crimped state when λ ≤ λ w . Yet, by assuming that a crimped fibre is energy-and stress-free, the expressions for P r and L r reduce to the respective integral terms, because W f (λ; λ) = ∂ λ W f (λ; λ) = 0.

Suitable Waviness Distributions Functions: Bounds and Normalisation
This subsection discusses suitable bounds and normalisation conditions for the waviness distributions functions p w .
Equation (9) postulates the existence of a lower bound λ w for the waviness λ w . A physically meaningful requirement for the lower bound is λ w > 0, since stretches are always nonnegative. Interestingly, this constraint allows the occurrence of fibres that are uncrimped and store elastic energy in the reference configuration. Thus, employing distribution functions p λw with lower bound λ w ∈ (0, 1) leads to the theory of residually stressed (or prestressed) materials, see, for example, Hoger [59] and Ogden [98] for details on these theories. In turn, requiring an energy-free (and stress-free) reference state of W r leads to the constraint λ w ≥ 1 for the lower bound, meaning that all fibres are crimped (or at least not actually stretched) in the reference state. There is no similar constraint for the upper bound λ w . However, we postulate that the waviness of the curliest fibre is still finite, since the only scenario in which this is not the case is when the two end-points of a fibre coincide and l 0 becomes zero, see Eq. (8). Hence, the set of continuous, univariate probability distributions with support R w = [λ w , λ w ] is appropriate for the waviness distribution function p w . Therein, the support defines the range in which the PDF is non-zero and positive, so that R w = {λ ∈ R : p w (λ) > 0}, whereas the PDF is zero if λ / ∈ R w . One prominent example for such a PDF is the four-parameter Beta distribution that has been used in the studies of, for example, Sverdlik and Lanir [118], Lokshin and Lanir [81], Chen et al. [27], Fan and Sacks [36], Fata et al. [37], Weisbecker et al. [127], and Avazmohammadi et al. [6]. It is introduced in detail in Appendix B, along with three alternative probabilities that have support on a bounded interval. However, since the definition of the upper bound is not as physically strict as the lower bound, also probability distributions with support on a semi-infinite interval [λ w , ∞) can prove useful and have been used in the literature, for example, the Weibull distribution by Hurschler et al. [63], the Gamma distribution by Sacks [109] and Bischoff [14] and the log-logistic distribution by Zulliger et al. [131] and Roy et al. [107]. The early works of Lake and Armeniades [71], Decraemer et al. [29], Lanir [73], Horowitz et al. [61] and Belkoff and Haut [12] considered a classical Gaussian normal distribution, which however may include unphysical non-zero probabilities for negative threshold stretches λ w due to its unbounded support and should thus be avoided.
In order to formulate a normalisation condition for p w , it is sensible to introduce the cumulative distribution function (CDF) The usual property of the CDF, F w (λ ≥ λ w ) = 1 for distributions on a bounded interval [λ w , λ w ], directly leads to the normalisation condition for p w : Corresponding conditions to Eqs. (14) and (15) for distributions on the semi-infinite interval [λ w , ∞) are defined via the limits lim λ→∞ F w (λ) = 1 and lim λ→∞ λ λ w p w dλ w = 1.

Calibration to Experimental Data
The choice of a specific waviness distribution function and the determination of the associated parameters requires suitable experimental data of the tissue that is modelled. One has to distinguish-roughly speaking-between two approaches: (i) Direct calibration of p w or F w to experimentally measured distribution data, such as histograms, or (ii) inverse estimation of p w or F w using statistical data in form of a discrete number of moments, such as the mean.   [104]. Those studies provide discrete, yet sufficiently high resolution data of the waviness (or straightness) distribution. Mohammadkhah et al. [91] provided waviness data in form of the moments (the mean and the standard deviation) that can be used for the second approach. Also note that this approach is reminiscent to the classical moment problem where the distribution must be determined on the basis of a discrete (and usually small) number of moments.
In contrast to these two bottom-up approaches, an alternative way is to choose an ODF in an ad-hoc manner and to determine the corresponding parameters by calibrating the resulting model to the overall fibre population or network behaviour, hence, calibrating P r orP , respectively. This has been done, for example, by Sverdlik and Lanir [118], Sacks [109], Weisbecker et al. [127] and Avazmohammadi et al. [6]. However, as the overall models usually result in a large number of parameters, it may be difficult with this method to establish a direct and unique link between the parameters obtained through calibration and the microstructural properties.
A direct calibration of the Beta, the Kumaraswamy, the triangular and the PERT distribution to the experimental waviness measurements of Roy et al. [107] and Chen et al. [26] is shown in Figs. 5 and 6, respectively. The experimental data and the calibrated PDF p w are shown on the left-hand side and the corresponding CDF F w on the right-hand side. The calibration was performed by means of the nonlinear least-squares method using the commercial software Matlab (https://www.mathworks.com). Note that the study of Roy et al. [107] provides the distribution in terms of the straightness, which however can be directly converted to waviness data. The three characteristic statistical moments, the mean, the standard deviation and the skewness, resulting from the best-fits in Figs. 5 and 6 are provided in Table 1. Figures 5 and 6 depict the asymmetric and right-tailed shape of the waviness distributions, resulting in a positive skew s w as noted in Table 1. This shape has also been observed in the experiments of Chen et al. [28] and Rezakhaniha et al. [104]. However, though very  This conclusion is drawn from taking the differences in the meanλ w and the standard deviation σ w between the measurements of Roy et al. [107] and Chen et al. [26]. Such structural heterogeneities among different specimen are typical for biological tissues, in which form and structure are highly optimised to the certain function of the material. Further, it can be observed that the Beta as well the Kumaraswamy distributions provide excellent fits to the experimental measurements and that they are very similar in their abilities and results. The triangular distribution is not able to match the experimental results as good as the other distributions, but still obtains for F w a curve that is very similar. Thus, it provides a very good balance between simplicity and accuracy. Table 1 underpins this statement by means of the statistical moments, as the best fit of the triangular distribution entails mean values almost equal to the more sophisticated Beta and Kumaraswamy distributions and is further able to predict the positive skew.

Modelling of Individual Crimped Fibres
The structural modelling approach for fibre populations relies on the statistical description of the waviness by means of the distribution function p w and on the mechanical description of single fibres by means of a strain-energy function W f (λ; λ w ). Further, it has been stated in Sect. 3.1 that the fibre energy in a stretch-based framework depends on the stretch λ and contains the waviness λ w as a parameter. This is essential to distinguish between the different fibre deformation states that were presented in Fig. 4. Particularly, fibres-such as collagen or elastin fibres in soft tissues-usually show a resistance to stretches λ > λ w that is magnitudes higher than in the crimped state when λ < λ w . A useful approach to account for this fact is to split up the energy W f into two terms via Therein, W cb f is the energy during compression and (un)bending, whereas W s f is the energy when the fibre is stretched beyond the waviness λ w . Models that describe the bending stiffness in the crimped state (hence W cb f ) were formulated by, for example, Freed and Doehring [42], Garikipati et al. [46], Grytz and Meschke [52] and Marino and Wriggers [83]. However, fibres in soft tissues are commonly idealised as structures without compressive and bending stiffness, see, for example, Lanir [72] and Kastelic et al. [67]. This assumption entails that the fibres do not store elastic energy and are stress-free whenever λ < λ w such that W cb f = 0. This leads to a λ-P f -diagram as depicted in Fig. 3. In this case, the fibre strain energy is fully defined in terms of the stretch-based contribution such that W f = W s f . Further, appropriate strain measures are needed that quantify the stretch after uncrimping 1 in order to formulate strain-energy functions that are zero whenever λ < λ w . To do so, a useful concept is the introduction of a multiplicative decomposition of the deformation gradient via where F w : r 0 → r w = F w r 0 describes the deformation from the reference configuration to the stress-and energy-free intermediate configuration. The decomposition concept was originally introduced by Lee and Liu [79] as F = F e F p in the context of elastoplasticity. However, in contrast to the original usage where F p and the corresponding intermediate configuration described the plastic (or inelastic) deformation that is determined by the motion history of the body, the tensor F w in Eq. (17) is different in nature. The intermediate configuration according to (17) is defined such that the stretch λ of a single fibre is equal to the waviness, hence, λ w = |r w | = √ F w r 0 · F w r 0 and F w therefore has the characteristic of an a priori known material property. However, the decomposition in (17) is not unique, since it is only defined by the one-dimensional stretch λ w . Any transformation QF w in which Q ∈ MG rw ti belongs to the transversely isotropic symmetry group MG rw ti (e.g., with arbitrary rotations around r w ) does not influence the results. Furthermore, the elastic stretch, λ e , of a single fibre after uncrimping is defined as the ratio between the length |r| = √ F r 0 · F r 0 = λ in the current configuration with respect to the length |r w | = λ w in the intermediate configuration, hence, λ e = λ/λ w . Consequently, we obtain as a one-dimensional version of (17) mapped in fibre direction. Further note that the elastic deformation part, F e : r w → r = F e r w , maps the intermediate vector r w to the current configuration, so that |F e r w | = λ.
It is now possible to define different strain measures to quantify the stretch deformation after uncrimping. Here, the aim is to introduce appropriate one-dimensional versions of the Lagrangean strain tensor family which is commonly referred to as Seth-Hill family due to Seth [112] and Hill [58] (but the concept has already been introduced by Doyle and Ericksen [30]). In Eq. (19), λ i are the principal stretches, n (i) are the corresponding eigenvectors and i = 1, 2, 3 in the present three-dimensional setup. Further, the principal stretches are the eigenvalues of the right stretch tensor 2 U = √ F T F . The tensor family in (19) resembles the logarithmic Hencky strain for m = 0, the Biot strain for m = 1 and the Green-Lagrange strain for m = 2. For the case m = 0, the logarithm arises as the consistent limit for m → 0, hence, lim m→0 (λ m − 1)/m = ln[λ]. Now, different possibilities exist to define strain measures that account for the energy-and stress-free nature of the intermediate configuration. The different approaches have also been investigated in the context of elastoplasticity, see, for example, Simo and Ortiz [115], Miehe [86] or Miehe et al. [87]. In particular, two strain families are usually considered:

Material Strain Description
A family of elastic strain tensors which has been advocated by Green and Naghdi [50,51] in the context of elastoplasticity proceeds from the additive decomposition of the strain tensor into an elastic part E (m) e and an inelastic (here "crimped") part  For m = 2, the resulting elastic Green-Lagrange-type strain E (2) e has been proposed frequently in the literature and has been denoted, for example, as e KL by Green and Naghdi [50], as E by Green and Naghdi [51], as E (e) by Kleiber [68] and as E e by Simo and Ortiz [115]. The designation as material strain is motivated through the fact that the right stretch tensor is defined in the reference configuration and is thus a material quantity.
In order to obtain appropriate one-dimensional strain measures from Eq. (20), we replace the occurring stretch tensors by the respective maps in fibre direction such that U → |U r 0 | = |F r 0 | = λ and U w → |U w r 0 | = |F w r 0 | = λ w . This leads to a class of onedimensional strains in the form of where the expression for E (0) e is obtained by means of the standard calculation rule ln[λ] − ln[λ w ] = ln[λ/λ w ] for λ, λ w > 0. It is obvious that the strains E (m) e are zero whenever λ = λ w . Further, the derivative of the strains with respect to the stretch λ is given by Hencky-, Biot-and Green-Lagrange-type strain measures are obtained from (21) by setting m = 0, 1, 2, respectively, and read They are visualised in Fig. 7a for two different waviness parameters λ w = {1.0, 2.0}. For later use, the strain measures from Eq. (21) evaluated for λ w = 1 (thus, λ e = λ) are denoted as Hence, these are the one-dimensional versions of the strain tensors from Eq. (19).

Intermediate Strain Description
As noted before, the strain measures defined in Eq. (20) This motivates the designation as multiplicative strain approach. The strain tensor group is then given bỹ Therein, λ m e(i) are the elastic principal stretches (the eigenvalues of U e ) and n e(i) are the corresponding eigenvectors that are based in the intermediate configuration. For m = 2, the strain tensorẼ (2) e corresponds to the frequently used intermediate Green-Lagrange strain that was denoted, for example, by Lee and Liu [79] as E e , by Simo and Ortiz [115] asĒ e and by Haupt [54] asˆ e andÊ e .
The respective one-dimensional strain measures are derived by replacing in Eq. (25) the elastic stretch tensor and the identity tensor with their one-dimensional counterpart, hence, U e → λ e and I → 1. This gives the strain measures The derivative with respect to the stretch λ is given by Hencky-, Biot-and Green-Lagrange-type strain measures are obtained from (26) by setting m = 0, 1, 2, respectively, and read The latter strain,Ẽ (2) e , has been denoted by Lanir [74] and Lanir [76] as e t and e f , respectively. The three strains from Eq. (28)

Fibre Energies with Proportional Stress-Strain Relationship
The formulation of strain-energy functions for initially crimped fibres includes two major steps: (i) The calibration of a suitable energy function (a mathematical form) to experimental data of uncrimped fibres that are subjected to uniaxial stretch and (ii) the choice of an appropriate strain description, like one of the previously presented ones, in order to account for the initially crimped state of the fibre and to include the waviness as a parameter. The first step is thereby motivated through the fact that in an experimental setup, single fibres are usually pre-stretched until they are uncrimped and straight. Experimental measurements proceed from the uncrimped state, hence λ w = 1. This means that a strain-energy function is formulated by means of "ordinary" strain measures E (m) and calibrated to experimental data of initially uncrimped fibres. Subsequently, the strain measure E (m) is replaced by the respective (with same order m) additive strain E (m) e or multiplicative strainẼ (m) e . Hence, the calibration step leads to the same material parameters, regardless of the choice of additive or multiplicative strains, but the strain energy will be different for waviness values λ w > 1. However, it has already been noted that for m = 0, the resulting logarithmic strain measures, E (0) e andẼ (0) e , are identical and make the choice between additive and multiplicative approach superfluous. Due to this interesting feature, we will provide more detailed investigations on logarithm-based energies in the next section.
Here, we focus on strain energies for single fibres that contain one stiffness parameter, μ f , and show a direct (linear) proportionality between strain and stress. This is motivated by the fact that fibres in soft biological tissues often exhibit an almost linear relationship between applied strain and the resulting stress. However, while this proportionality is unique in linear elasticity, it is not true for the finite-elasticity case. In particular, there exists no Young's modulus as proportionality constant at finite deformations, because it has to be specified which strain and which stress are linearly connected. Only experimental investigations can clarify this ultimately. Specific examples and calibrations to experiments are shown in Sect. 6, while the focus here is on a generic introduction and comparison of appropriate energy formulations.
In finite elasticity, it is sensible to establish a direct proportionality between, for instance, the experimentally measurable first Piola-Kirchhoff or the Cauchy stresses on one side and Hencky, Biot or Green-Lagrange strains on the other side. While the first Piola-Kirchhoff stress has already been introduced in Eqs. (10) and (11) and relates the current force to the referential (undeformed) geometry, the Cauchy stress (often referred to as true stress) relates the current force to the current (deformed) geometry. The scalar-valued first Piola-Kirchhoff stress of a single fibre is given as and enters the fibre population stress via Eq. (11), whereas the scalar Cauchy fibre stress is given by σ f = λ P f . A frequently used fibre formulation establishes a direct proportionality between first Piola-Kirchhoff stress and Biot strain via meaning that P P -E (1) f ∝ μ f E (1) . Therein, ∝ μ f denotes a direct proportionality with constant μ f . Equation (31) induces the quadratic strain energy that has been used, for example, by Humphrey et al. [62], Alastrué et al. [1] and Ghaemi et al. [48]. Note that as an alternative to the case distinction in the energy formulation, the distinction could also be applied directly to the strain measures (such that E (1) = 0 if λ < 1) as it was done by, for instance, Lanir [76]. At this point, two versions of (32) can be formulated for initially crimped fibres, that is, replacing E (1) by E (1) e orẼ (1) e . The first choice leads to the energy and the first Piola-Kirchhoff stress This additive formulation has been used by, for example, Lake and Armeniades [71], Cacho et al. [21], Chen et al. [27] (with an additional cubic energy term), Weisbecker et al. [127] and Bleiler et al. [16]. In contrast, the second, multiplicative strain leads to the energy and the first Piola-Kirchhoff stress It is interesting to note that while the additive strain formulation P e , this is not the case for the multiplicative formulation P P λw-Ẽ (1) e f , for which the proportionality constant is μ f /λ w , such that e /λ w . For further investigations, it is useful to consider the second derivatives of the energy functions, which denote the nominal stiffnesses. They are given as and  Figure 8a depicts that the stiffness based on the additive formulation exhibits constant values independent of the fibre waviness, while the multiplicative formulation shown in Fig. 8b proceeds from different stiffness values depending on the initial fibre crimp. More specifically, the stiffness of a single fibre decreases with increasing waviness and remains constant after recruitment. Now, as the definition of stiffness moduli in finite elasticity is not unique, a so-called incremental modulus L f can be formulated in addition to the nominal modulus L f , see, for example, Mihai and Goriely [89]. It is defined as and depicts the derivative of the Cauchy (true) stress with respect to the logarithmic (true) strain. 3 For the two strain energies defined in Eqs. (33) and (35), the incremental moduli are given by and The two moduli are plotted in Fig. 9 for different waviness values λ w = {1.0, 1.2, . . . , 2.0} and normalised by μ f , hence, for the same scenario as shown in Fig. 8. It can be observed that such that the first Piola-Kirchhoff stress is given by This stress formulation has been used by, for example, Decraemer et al. [29], Lanir [73,75] and Hurschler et al. [63]-without explicit formulation of the energy in (42). The corresponding stiffness moduli are given by as well as At this point, we note again that the three strain energies in Eqs. (33), (35) and (42) all reduce to Eq. (32) for λ w = 1. Hence, a calibration to experimental data of initially uncrimped fibres leads to the same stiffness parameter, μ f , for each energy, but the different energy formulations subsequently entail different behaviours for waviness values λ w > 1.
Alternatively to the first Piola-Kirchhoff stress, a direct proportionality can be established between a specific strain and the Cauchy stress. Doing so, one can formulate the energy such that the corresponding Cauchy stress is directly proportional to the Biot-type multiplicative strainẼ (1) w : and σ e . The stress formulation given in Eq. (47) has been proposed by Lanir [76]. Unlike the so far considered energies, the one presented in Eq. (46) is not quadratic. Further, the formulation is identical to the energy function proposed by Markert et al. [84] (when choosing therein M f = 1 and γ 1 = 1), which was inspired by the polynomial energy forms of Ogden [95]. The nominal stiffness due to the energy in (46) is given by and the incremental modulus is obtained as A visualisation of the two moduli for different waviness parameters is provided in Fig. 11. It is interesting to note from Fig. 11a that the nominal stiffness L σ -Ẽ (1) e f shows strain weakening and that the moduli of fibres with different waviness are on the same line. Further, as can be seen in Fig. 11b, the incremental stiffness L σ -Ẽ   Table 2 summarises the characteristic stress and stiffness values of the four energy functions presented in Eqs. (33), (35), (42) and (46). Many more energy functions with linear stress-strain relationship can be formulated in the framework of finite elasticity, for example, an energy W This formulation produces a direct proportionality between the second Piola-Kirchhoff stress, defined as S f = ∂ E (2) W f = λ −1 ∂ λ W f , and the Green-Lagrange-type multiplicative strainẼ (2) e and has been used, for example, by Sverdlik and Lanir [118], Lokshin and Lanir [81] (together with an additional nonlinear term) and Sacks [109]. Fan and Sacks [36] and Fata et al. [37] formulated the similar quadratic energy W Sλ 2 w -Ẽ (2) e f = 1 2 μ f (Ẽ (2) e ) 2 that entails the waviness-dependent proportionality constant μ f /λ 2 w between the second Piola-Kirchhoff stress and the multiplicative Green-Lagrange-type strainẼ (2) e . Again, these two energies, (2) e f , are identical for λ w = 1 and yield the same fibre stiffness μ f if calibrated to data of uncrimped fibres, but they will be different whenever λ w > 1. Finally, the Table 2 Stresses and stiffness quantities for four exemplary fibre energies (for λ ≥ λ w ) Biot-type multiplicative strainẼ (1) e was applied by Bircher et al. [13] and Li and Holzapfel [80] in order to incorporate the description of initially crimped fibres into the energy formulations of Rubin and Bodner [108] and Holzapfel et al. [60], respectively.

Energies Based on Logarithmic Strains
As already mentioned in the previous section, the Hencky-type logarithmic strain ln[λ e ] = E (0) e =Ẽ (0) e is a very appropriate measure for the quantification of the stretch after uncrimping, because it makes the choice between additive (E (0) e ) and multiplicative (Ẽ (0) e ) strain formulation superfluous. The advantage of Hencky-type strains in connection with a multiplicative decomposition of the deformation, as formulated in Eq. (17), is already known for a long time. For instance, Bazant [9] formulated the key feature that "after one deformation [here the stretch λ w ], the new configuration can be taken as the reference state for computing the additional strain for a further deformation". Neff and Ghiba [92] gave similar explanations about the advantages of logarithmic strains in the context of elastoplasticity. A quadratic potential in the logarithmic Hencky-type strain ln[λ e ] = E (0) e =Ẽ (0) e reads W H f = μ f (ln[λ e ]) 2 /2 and represents a one-dimensional version of the classical Hencky strain energy-supplemented by the ability to describe initially crimped fibres. The corresponding Cauchy stress σ H f = μ f ln[λ e ] is directly proportional to the logarithmic strain, with σ H f ∝ μ f ln[λ e ], and implies the interesting property that the corresponding incremental modulus L H f = μ f becomes a stretch-independent constant. The Hencky energy is therefore considered as the energy that resembles linear elasticity in the finite deformation regime very closely. However, it is well known that the classical Hencky energy violates the convexity requirement at large strains, see, for example, Neff et al. [93,94]. To circumvent this shortcoming, Neff et al. [93] and Schröder et al. [111] introduced the so-called exponentiated Hencky energy, followed by Govindjee et al. [49] who formulated a corresponding one-dimensional version that can be used as fibre energy. A version that incorporates the initial crimp of single fibres consequently reads where k is dimensionless parameter that governs the nonlinearity. The corresponding first Piola-Kirchhoff stress is then given by with a nominal modulus and an incremental modulus The two moduli are visualised in Fig. 12 for k = 1 and λ w = {1.0, 1.2, . . . , 2.0}. The nominal moduli shown in Fig. 12a are non-monotone, unlike the constant or monotone behaviours of the functions investigated so far. This means that an initial strain softening is followed by a strain stiffening at higher strains. In contrast to that, the incremental modulus in Fig. 12b shows a monotone behaviour and high strain stiffening at higher strains. However, there is no initial stiffening of the incremental modulus, meaning that ∂ λ L eH f | λ=λw = 0. This resembles the linear characteristic of the classical Hencky energy.
At this point, we remark that the usage of Hencky-type energies usually requires the (manageable, but yet slightly cumbersome) derivation of eigenvectors for the computation of stress or elasticity tensors, see, for example, Ogden [97] or Simo [114]. This does not apply for the here investigated one-dimensional fibre energy, since the strain ln[λ e ] = E (0) e =Ẽ (0) e is a direct function of λ (and thus of F ) and not of ln [U ]. Hence, no polar decomposition is required.
In addition, a strain tensor family that approximates the logarithmic strains was proposed by Bazant [9].
respectively. This allows the formulation of a variety of alternative energy functions with specific linear stress-strain relations and again illustrates that there is no unique Young's modulus in finite elasticity to describe the direct proportionality between stress and strain.

Novel Analytical Expressions for Fibre Populations
The structural modelling of fibre populations presented in Sect. 3.1 results in integral expressions for the energy W r in Eq. (9) and, consequently, for the stresses and elasticity moduli based thereon. In turn, these effective fibre population quantities enter the effective network quantities that are integral expressions as well. For instance, the effective network energȳ W V in terms of a Voigt-type averaging is derived from Eq. (7). The integral over the microsphere surface, as defined in Eq. (4), can usually not be solved analytically and numerical quadrature schemes have to be applied. This means that a sufficient number of quadrature points is required to guarantee the accuracy of the schemes, see, for example, Ehret et al. [33], Verron [123] or Itskov [64]. Let N be the number of quadrature points, then the integral for W r has to be evaluated N -times to obtainW V . Thus, analytical expressions for the fibre population energy W r and its derivatives potentially lead to a significant reduction in computational time. While solving the integrals in Eqs. (9), (11) and (13) generally has to be done in a numerical manner, it is possible to formulate analytical expressions for specific choices of W f and p w . This is done in this section for fibre energies based on the additive Biot-type strain E (1) e . Such energies allow to write the integrals in Eqs. (9), (11) and (13) as convolutions, which proves useful for finding analytical solutions by means of symbolic computation software like Mathematica (https://www.wolfram.com/mathematica/).

Convolution Integrals for Effective Fibre Population Quantities
If the fibre energy W f (E (1) e ) depends on the additive Biot-type strain E (1) e = λ − λ w , then the effective fibre population energy in Eq. (9) is written in terms of a convolution by Strictly speaking, the convolution is equivalent to the integral with integration boundaries from −∞ to λ. Yet, since the value of the integral is zero from −∞ to λ w , the convolution is equivalent to the integral. Subsequently, with the derivation rule ∂ x (f * g) = (∂ x f ) * g = f * (∂ x g) for two arbitrary functions f and g, the integral formulations for the stress and the stiffness in Eqs. (11) and (13) can be written as with P f = ∂ λ W f , see Eq. (30), and L f = ∂ 2 λ 2 W f . In the remainder of this section, analytical expressions are presented for fibre energies of the frequently used form presented in Eq. (33), hence, W . It has been shown that the Fig. 13 The sharp stiffness jump of a single fibre at straightening and the effective nominal stiffness of a fibre population with beta-distributed waviness, where the variance σ 2 w has the numerical character as a regularisation parameter thereon based nominal stiffness, given in Eq. (37), denotes a step function with step μ f at λ w . With this, Eq. (55) 2 simplifies to The CDF, F w , of waviness distribution functions is typically of sigmoidal shape and is used as regularisation of sharp jumps also in other applications. This shows that beyond the physical meaningfulness, the structure of fibre population energies, together with fibre energies based on the additive strain E (1) e , smooths the sharp stiffness jump of single fibres at straightening and serves as a sensible mathematical and numerical regularisation. This is illustrated in Fig. 13.

Novel Analytical Expressions for Beta-Distributed Fibre Waviness
If single fibres are described by Eq. (33), hence, W P -E (1) e f , and the waviness of the fibres is beta-distributed, such that λ w ∼ Beta(β 1 , β 2 , λ w , λ w ), the effective energy of a fibre population is given in the analytical form where Therein, 2 F 1 is the Gaussian hypergeometric function, 2F1 is the regularised Gaussian hypergeometric function, [c] = (c − 1)! is the Gamma function and B is the Beta function.
More details on such mathematical functions can be found, for example, in the mathematical handbooks of Bronshtein et al. [19] and Olver et al. [99]. The Gaussian hypergeometric as well as other herein occurring functions are implemented in mathematical software tools like Matlab or part of numerical libraries like SciPy (https://www.scipy.org). They are thus easily accessible. Note, the energy in Eq. (57) has been presented in a slightly different notation by Cacho et al. [21]. Subsequently, the scalar stress P r is obtained either by means of a derivation of Eq. (57) with respect to stretch λ or by performing the convolution (55) 1

. In any case, it results in
and is the incomplete regularised Beta function that, in turn, contains the incomplete Beta function B x , see also Eq. (79). Further, the nominal stiffness is obtained by deriving Eq. (59) with respect to λ or by making use of the simple expression (56), yielding

Novel Analytical Expressions for Kumaraswamy-Distributed Fibre Waviness
If single fibres are described by Eq. (33) and the fibre waviness follows a Kumaraswamy distribution, then the effective energy of a fibre population is given by The herein occuring mathematical functions and expressions have been explained in the previous section.
The associated stress is given by and the nominal stiffness by

Novel Analytical Expressions for Triangularly-Distributed Fibre Waviness
For fibres that are described by Eq. (33) and exhibit triangularly-distributed waviness, λ w ∼ Tri(λ w , λ m w , λ w ), the effective energy is where the meanλ w and the variance σ 2 λw are given in Eq. (87). The associated stress is given by (67) and the nominal stiffness is given by This section is concluded by two remarks.
Firstly, we note that the analytical formulations for energy, stress and stiffness presented in this section are based on a case distinction to identify the different stretch regions. Similar energy forms based on such a case distinction were presented, for example, by Puso and Weiss [103] and Röhrle and Pullan [106]. However, the results obtained here are fully based on principles of structural modelling and not in an ad-hoc manner. The formulations directly result from the properties of single fibres and their structure by means of the waviness distribution. Moreover, the expressions for the stress P r resemble a bilinear stress-relation where the first linear region (if λ < λ w ) has negligible stiffness (since the energy W cb f is taken as zero) and the final stiffness (if λ > λ w ) is μ f . If the formulations are supplemented by an additional energy term, such as, for example, μ g (λ − 1) 2 /2 with μ g < μ f , then one obtains the bilinear forms that were presented by, for example, Babu et al. [7] and Freed and Rajagopal [43], here with first linear stiffness μ g and second linear stiffness μ g + μ f .
Secondly, we point out that the herein presented analytical results for the integral expressions at the fibre population level provide a significant reduction of computational costs compared to numerical solution schemes. We measured that the difference between the evaluation of the analytical results and adaptive as well as trapezoidal integration schemes is at least one magnitude. At the fibre network level, where effective population values have to be computed multiple times, we noticed a performance advantage of at least two magnitudes. For the sake of brevity, the associated details and results can be found in Appendix C.

Examples and Discussion
Next, we present examples and case studies in order to compare the formulations presented so far. Particularly, we start in Sect. 6.1 with a calibration of the different energy formulations for single collagen fibres from Sects. 4.3 and 4.4 to experimental data. Those results are employed in Sect. 6.2 to compare the different fibre energy formulations that were presented in Sects. 4.3 and 4.4 and hence the different strain measures that were presented in Sects. 4.1 and 4.2 on the level of fibre populations. Subsequently, the different fibre energy formulations are compared on the level of the entire fibre network. This is done by making use of the Voigt-based formulation that was formulated in Sect. 2.2.
In order to draw meaningful conclusions, we will only consider physiologically sensible scenarios in this section. We will therefore not choose generically selected parameters for which the differences between the formulations can be made arbitrarily large in an artificial way.

Calibration to Collagen Fibril/Fibre Experimental Data
In this subsection, experimental data on single, initially uncrimped collagen fibrils/fibres is employed to calibrate three of the strain-energy functions from Sects. 4.3 and 4.4. In this regard, it should be noted that collagen fibres are bundles of densely packed collagen fibrils, see, for example, Fratzl [41]. In turn, collagen fibrils are themselves fibre-like structures with a diameter in the range of 50-500 nm [see, e.g., 100,128]. Likewise, the diameters of collagen fibres that were reported in the literature range between 5 µm up to 500 µm [see, e.g., 3,35,47] and thus already reach the magnitudes of, for example, skeletal muscle fibres. The separation between fibrils and fibres in mechanical models is not always sharp and the term "collagen fibre" might generally refer to a fibrous collageneous structure in soft tissues. Numerous studies have dealt with the experimental investigation of the stressstrain behaviour of collagen at different scales, see, for example, Sasaki and Odajima [110], Fig. 14 Best-fit stress curves of three fibre energy formulations, calibrated to the experimental collagen fibril data of Sasaki and Odajima [110] Miyazaki and Hayashi [90], Gentleman et al. [47], Bozec and Horton [18], van der Rijt et al.
We use the energy W P -E (1) f from Eq. (31) that resembles a linear relationship between first Piola-Kirchhoff stress and Biot strain. Note that for the considered scenario of initially uncrimped fibres (with λ w = 1), the energies W , as defined in Eqs. (33), (35) and (42), respectively, are equal to W P -E (1) f . Further, the energy W σ -Ẽ (1) e f from Eq. (46), which constitutes a linear relationship between Cauchy stress and Biot strain and the exponentiated Hencky energy W eH f , is considered (with λ w = 1). The results of a least-squares calibration of the first Piola-Kirchhoff stresses, derived from the three considered energies, to the experimental data of Sasaki and Odajima [110] are presented in Fig. 14 the coefficients μ f = 441.4 MPa and k = 53.1 with R 2 = 0.9492 (however, the confidence interval for k is very large and choosing k = 1 results in a fit with R 2 = 0.9477 and is thus virtually equally well). It can be seen from Fig. 14 that the experiments of Sasaki and Odajima [110] revealed a linear relationship between the first Piola-Kirchhoff stress and the applied stretch. However, in the quite small stretch range that was investigated, the three energies result in very similar stiffness parameters μ f and are hardly distinguishable. Concluding, each of the considered energy functions is suitable for the description of single collagen fibrils.

Comparison of Fibre Energies and Strain Measures at the Fibre Population Level
Effective quantities of fibre populations are considered in this subsection in order to compare the different strain measures and fibre energies from Sect. 4. For the sake of clarity, we focus on two scenarios. The first one considers the fibre population with beta-distributed waviness λ w ∼ Beta(1.398, 3.523, 1.013, 1.099) as measured by Roy et al. [107], see Fig. 5 and the second such with a beta-distributed waviness λ w ∼ Beta(3.572, 17.89, 1.016, 2.077) measured by Chen et al. [26], see 6. Hence, the first scenario considers that the distribution of the fibre waviness is rather in the small-stretch regime while the second one leads to a wider distribution also at higher stretches. Further, we employ five different energy formulations for the single fibres, namely, W and W eH f , defined in Fig. 15 (a) Strain energy W r and (b) nominal stress P r of a fibre population with five different fibre energy formulations. The fibre waviness is beta-distributed according to experimental data of Roy et al. [107] and the fibre material parameters are chosen according to a calibration to the data of Sasaki and Odajima [110] Eqs. (33), (35) and (42), (46) and (50), respectively. We choose the material parameters that were obtained from the calibration to the experimental data of Sasaki and Odajima [110] in the previous subsection. Hence, we take μ f = 443.2 MPa for W f and μ f = 441.4 MPa for W eH f . Moreover, for the sake of comparability, we choose k = 1 for the exponentiated Hencky energy.
The fibre population energy W r and the corresponding nominal stress measure P r for the first scenario are shown in Fig. 15. The integrals in Eqs. (9) and (10)  (1) e f are nearly identical. Hence, in this scenario, the choice between additive and multiplicative strain measure is almost equivalent for the considered fibre strain energies that are based on first Piola-Kirchhoff stress and Biot strain.
This advocates the use of the additive-strain-based formulation, W P -E (1) e f , for which analytical expressions were presented in Sect. 5. It further has to be emphasised that the difference between the three first-Piola-Kirchhoff-stress-based formulations and the remaining two energies does not allow for a judgement in the sense of accuracy, since the difference is simply based on different assumptions how the fibre behaviour proceeds beyond the stretch range for which they were calibrated. These different behaviours can also be seen in Fig. 16, in which the nominal stiffness measure L r and the incremental modulus L r are visualised. The first three formulations assume a constant nominal stiffness beyond the stretch at which all  The fibre waviness is beta-distributed according to experimental data of Chen et al. [26] and the fibre material parameters are chosen according to a calibration to the data of Sasaki and Odajima [110] fibres are recruited, whereas the Cauchy-stress-based formulation W σ -Ẽ (1) e f shows a nominal strain weakening and a linear strain stiffening in the incremental modulus. The exponential Hencky energy W eH f first shows moderate strain stiffening of the incremental modulus, see Fig. 16b, and shows increasing stiffening at higher stretches due to the nonlinear nature of the W eH f (even for k = 1). The energy and the nominal stress for second scenario, hence, with a wider fibre waviness distribution λ w ∼ Beta(3.572, 17.89, 1.016, 2.077) according the experiments of Chen et al. [26], are shown in Fig. 17. It can be seen that the results are similar to ones shown in Fig. 15 from a qualitative point of view. Yet, the three first-Piola-Kirchhoff-stress-based formulations differ more than in the previous scenario. This shows that the difference between additive and multiplicative strain formulation comes more into play for wider waviness distributions. However, the difference only becomes important at high stretches, whereas the stress responses based on the five fibre energy formulations are fairly similar up to a stretch

Comparison of Fibre Energies and Strain Measures at the Fibre Network Level
In the last example, we examine the different fibre energy formulations, W and W eH f , at the level of the entire fibre network. The material parameters for the fibre energies are chosen according to the calibration to the data of Sasaki and Odajima [110] as obtained in Sect. 6.1. To this end, we examine four different scenarios resulting from the combination of the two fibre waviness distributions from the previous subsection, λ w ∼ Beta(1.398, 3.523, 1.013, 1.099) and λ w ∼ Beta(3.572, 17.89, 1.016, 2.077), and two macroscopic deformations. The first deformation is an incompressible 4 uniaxial extension that is described by the macroscopic deformation gradientF = 1/ √ λ (e s 1 ⊗ e s 1 + e s 2 ⊗ e s 2 ) + λ e s 3 ⊗ e s 3 with λ ∈ [1.0, 2.0] and the second is a simple shear deformation described bȳ F = I + γ e s 3 ⊗ e s 1 with γ ∈ [0, 1.5], where I is the second-order identity tensor. We assume a uniform fibre distribution over the microsphere surface so that the ODF is given by p = 1/(4π). Further, for the sake of simplicity, we provide the effective network response on the basis of the affinity assumption as presented in Sect. 2.2. The surface averaging operator defined in Eq. (4) is calculated by means of the numerical quadrature scheme by Lebedev [77] and Lebedev and Laikov [78] and we choose N = 770 quadrature points to guarantee sufficient accuracy.
The first scenario considers a fibre waviness distribution λ w ∼ Beta(1.398, 3.523, 1.013, 1.099) due to the experiments of Roy et al. [107] and a uniaxial extension of the fibre network. The resulting effective (Voigt-type) energiesW V , defined in Eq. (7), and the thereon based nominal stresses ∂ λWV are visualised in Fig. 19. The curves of the effective energȳ W V in Fig. 19a resemble to the ones presented in Fig. 15a for the fibre population energy are significantly lower than the other three ones, which also been observed in Fig. 15a. Moreover, it can be seen from the magnified view in Fig. 19b that the effective responses, in this case the stress responses, are virtually identical up to a stretch λ = 1.1.
In the second scenario we investigate the same uniaxial deformation as in the first one but with the wider fibre waviness distribution λ w ∼ Beta(3.572, 17.89, 1.016, 2.077) based on the experimental data of Chen et al. [26]. The results for the effective energyW V and the nominal stress ∂ λWV are shown in Fig. 20. It can be seen that the wider waviness distribution causes larger differences between the additive-strain-based formulation W P -E   Figure 21 shows the results for the fibre waviness distribution λ w ∼ Beta(1.398, 3.523, 1.013, 1.099) and results for the wider waviness distribution λ w ∼ Beta(3.572, 17.89, 1.016, 2.077) are visualised in Fig. 22. The results of the simple shear deformation are essentially similar to those previously considered. Particularly in Fig. 21 the first three fibre energy formulations lead to effective energies that are virtually identical, and W eH f , for their part, lead to effective energies that are close to each other and significantly lower than the other three formulations.

Concluding Remarks
In this paper, we discussed the modelling of fibre networks with a special focus on the statistical description of populations of unidirectionally-oriented and initially crimped fibres. In addition to providing an overview of the basic model equations and principles in Sects. 2 and 3, we have dealt with existing gaps and aspects of the topic that have not been discussed in detail so far. Particularly, we (i) formulated on the basis of the multiplicative deformation decomposition two strain families to describe the stretch of a fibre after uncrimping, (ii) demonstrated the different implications on the resulting stiffness-such as waviness-dependent or waviness-independent-of commonly used fibre strain energies, (iii) formulated an exponentiated Hencky energy for crimped fibres and (iv) presented easy-touse analytical expressions for fibre populations with beta-distributed waviness that make the-potentially expensive-calculation of integrals superfluous.
This work and the herein presented examples provided, among theoretical insights, the following interesting conclusions: • The two very commonly used fibre energy functions W P -E Hencky energy W eH f are very suitable for the description of crimped fibres and make the choice between additive or multiplicative strain measures superfluous. However, it must be assessed whether the specific (and potentially very nonlinear) fibre behaviour described by W eH f accurately describes the real behaviour of fibres under large deformations.

Appendix A: Effective Network Response via Variational Energy Minimisation
In Sect. 2.2, the affine map (5) leads to a Voigt-type formulation for the effective network re-sponseW V . For completeness, we provide some details on alternative nonaffine approaches in this appendix. Hence, in line with classical homogenisation treatments, see Hill [57] or Ponte Castañeda and Suquet [101], the micro-deformations F r of space orientations in network models are such that Ω is minimised and the effective energyW is defined in terms of a minimum energy principle asW Hence, the effective energy is a function of a macroscopic deformation gradientF that acts on the microsphere. The micro-deformations are linked to the macroscopic deformation gradient viaF and it is useful to formulate the variational problem (69) in terms of the Lagrange functional where P r is a tensor-valued Lagrange multiplier that enforces condition (70), see also Govindjee et al. [49]. Standard variational calculus demands that the variation δL vanishes and that the Lagrange functional becomes stationary with respect to F r and the Lagrange multiplier P r . While the stationarity condition of L with respect to P r results in Eq. (70) (as expected), stationarity of L with respect to the micro-deformations F r delivers The last expression represents the effective first Piola-Kirchhoff stress, as it is the average of the microscopic first Piola-Kirchhoff stresses P r := ∂ Fr W r (F r ). As a matter of fact, it is easy to show that the Lagrange multiplier tensor P r represents the macroscopic stress tensor P = ∂FW . By utilising the derivative of Eq. (70) with respect toF , giving the fourth-order identity tensor I = ∂F F r Ω , and assuming that F * r is a solution for the variational problem, we obtain Network models with effective properties on the basis of a minimisation principle, as formulated in (69), were presented by Govindjee et al. [49] and Chen et al. [27]. The latter proceeded from a discrete number of space orientations instead of a continuous distribution and formulated an estimate in terms of the tangent-second-order homogenisation method by Ponte Castañeda and Tiberio [102] [see also 17]. Moreover, Miehe et al. [88] formulated the so-called non-affine microsphere model on the basis of a variational minimisation principle for the scalar fluctuations of the stretch λ r and the tube contraction ν r on the surface Ω. An energy relaxation based on a so-called maximal advance path constraint was formulated by Tkachuk and Linder [119].
function (PDF) and cumulative distribution function (CDF), also characteristic measures that describe the statistics of the distributions are provided. These are the first raw moment which is the mean value (the expectation value) of the distribution, the second central moment which is the variance of the distribution, and the third normalised moment which is the skewness. In Eq. (76), the standard deviation σ w is the square root of the variance σ 2 w .

B.1 Beta Distribution
The four-parameter Beta distribution is widely used for the description of collageneous tissues. References were given in Sect. 3.2. The four-parameter distribution is a shifted and scaled version of the ordinary two-parameter Beta distribution that has support on the range [0, 1]. A beta-distributed waviness is denoted as and thus depends on two shape parameters, β 1 and β 2 , and on the lower and upper bound. The ODF and the corresponding CDF read correspondingly. Therein, B[β 1 , β 2 ] is the Beta function and I τ λ (τ λ ) is the incomplete regularised Beta function that, in turn, contains the incomplete Beta function The incomplete regularised Beta function thereby entails the characteristic sigmoidal shape of the CDF. The meanλ w , variance σ 2 w and the skewness s w of the four-parameter Beta distribution are given byλ and respectively. The Beta distribution is symmetric whenever β 1 = β 2 , but also allows to describe asymmetric distribution shapes for β 1 = β 2 . If β 1 = β 2 = 1, it results in a uniform distribution of the values between λ w and λ w . Another special case is given for β 1 = β 2 = 1/2, resulting in a bimodal Arcsine distribution with modes at the boundaries of the support. Hence, the PDF p w obeys two maxima. Restricting to β 1 , β 2 > 1 leads to unimodal distributions with mode λ m w = λ w + (β 1 − 1)/(β 1 + β 2 − 2)(λ w − λ w ).

B.2 Kumaraswamy Distribution
The four-parameter distribution by Kumaraswamy [70] is less known but equally powerful as the frequently applied Beta distribution, see also Jones [66]. A Kumaraswamy-distributed waviness is denoted as and depends on two shape parameters, a and b, and the lower and upper bound. The PDF and the CDF are given by respectively. Note that τ λ has been defined in Eq. (79). The meanλ w , variance σ 2 w and the skewness s w of the four-parameter Kumaraswamy distribution are given bȳ respectively. Therein,λ w01 and σ w01 are the mean and the standard deviation, respectively, of a Kumaraswamy(a, b, 0, 1) distribution in the interval [0, 1]. Hence, It is interesting to note that while the Beta distribution contains the Beta function in the expressions for the PDF and CDF in Eqs. (78), but not in the moment-based measures in Eqs. (80), it is exactly the opposite for the Kumaraswamy distribution. Hence, the PDF and the CDF in Eqs. (82) only contain regular expressions while the moments in Eqs. (83) are dependent on the Beta function.
Like the Beta distribution, the Kumaraswamy distribution can account for very general asymmetric distribution shapes and contains the uniform distribution as special case when a = b = 1. It is unimodal when restricting to a, b > 1. However, equality of the shape parameters, a = b, does generally not entail symmetry of the distribution.

B.3 Triangular Distribution
The triangular distribution is very simple yet useful probability formulation. A triangularlydistributed waviness follows a Beta distribution λ w ∼ Beta(β 1 , β 2 , λ w , λ w ) where All others measures can be calculated from the corresponding expressions for the Beta distribution, such as, for example, the moments in Eqs. (80).

Appendix C: The Numerical Advantage of the Analytical Expressions from Sect. 5
This appendix provides a brief comparison between different solution methods for the integral expressions in Eqs. (9), (11) and (13), hence, Exemplary, we focus on the advantage of the analytical expressions (57), (59) and (62) so that p w describes a beta-distributed fibre waviness. To do so, we compare the computational times of two numerical integration schemes to the ones of the analytical expressions. Clearly, computational (wall-clock) times are highly dependent on hardware (platform), software (programming language) and implementation, so that it does not make sense to provide absolute values. Hence, the focus is on qualitative comparisons for the three different methods-the two numerical and the analytical-and we normalise the times of numerical methods to the analytical one. Consequently, numbers higher than 1 indicate higher computational effort compared to the analytical solution. The studies are done in Matlab code (version 2022a; www.mathworks.com). The times are measured with the built-in functions tic/toc. Numerical integration schemes usually proceed from a discretisation of the integration interval, [λ w , λ], and approximate the integral in Eq. (91) as V d r ≈ V r by means of the evaluation of the integrand at a finite number of points, say N . For our studies, we choose an adaptive integration with the built-in function integral and the trapezoidal rule by means of the built-in function trapz. Note, the analytical solution in Eq. (57) contains the hypergeometric function, 2 F 1 . To compute this function, we employed the function hypergeomq, 5 because the corresponding Matlab-built-in function hypergeom is not designed for performance. For the sake of brevity, we choose in the following a beta-distributed waviness λ w ∼ Beta(1.398, 3.523, 1.013, 1.099) as measured by Roy et al. [107], which is also used in the examples in Sect. 6. The fibre stiffness is μ f = 1.
To get an idea of the integration to be solved, Fig. 23a exemplary shows the integrand of Eq. (9), W f p w , for the stretch λ = 1.1. The red circles mark 200 discrete points that are chosen for the trapezoidal integration scheme in the following. This choice is motivated by the results of a convergence study that is presented in Fig. 23b. It shows the absolute value of the relative error, |(V d r − V r )/V r |, of the trapezoidal scheme with respect to the analytical (exact) solution in dependence of the number of integration points N . The convergence of energy (solid lines), stress (dashed lines) and stiffness (dotted lines) computation is depicted by means of a log-log plot for three different stretches, λ = {1.1, 1.5, 2.0} (blue, red and green curves, respectively). The figure shows that higher stretches demand for a higher number of  Runtime comparisons between analytical, adaptive (blue) and trapezoidal (red) solution strategies for the integrals in Eqs. (9) (dark colours), (11) (medium colours) and (13) (light colours) for three different stretches λ = {1.1, 1.5, 2.0}. The integrals are carried out 50,000 times. The times of the adaptive and trapezoidal schemes are normalised to the time of the analytical solutions so that 1 means equal performance and values larger than 1 indicate better performance of the analytical method integration points to maintain the accuracy. In the considered example, 200 points appear appropriate, providing an accuracy of around 10 −2 at moderate stretches and also guaranteeing an error below 10 −1 at high stretches. Note that lower errors demand for a considerable higher number of integration points, for instance, an error of 10 −6 requires around 10,000 points. To provide a fair comparison, we also choose 10 −2 for the tolerances (RelTol and AbsTol) of the adaptive integration scheme. Figure 24 presents the normalised computational times of adaptive (blue) and trapezoidal (red) schemes for energy (dark colours), stress (medium colours) and stiffness (light colours). Three different stretches, λ = {1.1, 1.5, 2.0}, are shown. In order to avoid the influence of individual runs, each of the bars is the averaged result of 50,000 runs. It can be  Fig. 24 that there is at least one magnitude difference between numerical and analytical methods and that the extra effort for the numerical schemes increases with increasing stretch. Further, it is evident that the difference in the stiffness computation is significantly larger and almost reaches two magnitudes. This is mainly attributed to the fact that the analytical expression (62) becomes very simple and does not require any evaluation for λ > λ w at all. We want to emphasise at this point that the chosen number of integration points for the trapezoidal scheme and the tolerances for the adaptive scheme are rather moderate. The computational benefit of the analytical expressions will further increase when more accuracy is required from the numerical methods.
To obtain effective values on the fibre network level, the effective values of a fibre population have to be computed multiple times. For instance, we applied in Sect. 6.3 a numerical quadrature scheme with N = 770 points to obtain the effective network energyW V . This means that the integral in Eq. (91) has to be evaluated N times to obtain W r at each quadrature point. To represent such a scenario, we finally compare the performance of the different solution strategies for the integrals in Eq. (91) for a series of stretch values. Hence, we compare the times needed to compute a stretch interval λ = [1.0, 1.1] with N λ = {100, 1000} points as well as the interval λ = [1.0, 2.0] with N λ = 100 points. The results are presented in Fig. 25. In order to avoid the influence of individual runs, each of the bars is the averaged result of 5,000 runs. It can be seen that the numerical schemes are at least two magnitudes (100 times) slower than the evaluation of the analytical results. The adaptive schemes is for some cases even three magnitudes slower. The computational advantage of the analytical expressions increases with increasing N λ and increasing maximum stretch. Note that the number N λ in this scenario is a representative of the number of quadrature points N . The schemes of Lebedev [77] and Lebedev and Laikov [78] (used in Sect. 6.3) provide up to N = 5810. Such a high number is required, for example, if the fibre network is highly anisotropic. Hence, such scenarios will further increase the benefit of the analytical results compared to numerical schemes. Code Availability Not applicable.

Conflict of interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.