A Tractable, Transferable, and Empirically Consistent Fibrous Biomaterial Model

Stochastic modeling is a useful approach for modeling fibrous materials that attempts to recreate fibrous materials’ structure using statistical data. However, several issues remain to be resolved in the stochastic modeling of fibrous materials—for example, estimating 3D fiber orientation distributions from 2D data, achieving the desired fiber tortuosity distributions, and dealing with fiber–fiber penetration. This work proposes innovative methods to (1) create a mapping from 2D fiber orientation data to 3D fiber orientation probability distributions, and vice versa; and (2) provide a means to select parameters de novo for random walks employing the popularized von Mises–Fisher distribution given that the desired tortuosity of the path is known. The proposed methods are incorporated alongside previously developed stochastic modeling techniques to simulate fiber network structures. First, fiber orientation distributions vary significantly depending on how a fibrous material is formed, and projection distortion affects the measurement of fiber orientation distributions when reported as 2D data such as histograms or polar plots. Relationships are developed to estimate 3D fiber orientation distributions from 2D data, accounting for projection distortion and the variety of orientation distributions observed in fibrous materials. We show that without correcting for projection distortion, fiber orientation distribution parameters could have errors of up to 100%. Second, in stochastic modeling, fiber tortuosity is usually treated with random walks, but no relationship is available for choosing random walk inputs to generate a desired fiber tortuosity. Relationships are also developed to relate the input parameters of von Mises–Fisher random walks to the expected tortuosity of the generated path—a necessary link to modeling fiber tortuosity distributions tractably and with empirical consistency. Using the developed relationships, we show that modeling of tortuous fibers from a distribution could be sped up by ~1200-fold and the uncertainty of selecting appropriate parameters could be eliminated. Third, randomly placing fibers in a simulation domain inevitably results in fiber–fiber penetration, and correcting this issue requires changes to the simulated fibrous material structure through non-penetration conditions. No thorough remedy can be offered here, but we statistically quantify the effects of enforcing non-penetration conditions on the fiber shape and orientation changes as well as the overall fibrous material model. This work offers tractable and transferable methods for treating fiber orientation and tortuosity that allow for empirical consistency in the stochastic modeling of fibrous materials.


Introduction
Numerical modeling is a valuable tool for investigating the properties of complex fibrous materials (FMs). For instance, numerical modelling has offered deeper knowledge regarding the mechanics of FMs embedded with stiff inclusions [1], the elasticity and plasticity of plant cell walls [2], the heat and momentum transfer between FMs and shear-thinning liquids [3], and the acoustic properties of FMs [4]. The drive to obtain more detailed insights and to understand the structure-property relationships of FMs naturally inspires researchers to develop more realistic FM models. Existing approaches to modeling FM structures include physics-driven approaches [5][6][7], random sequential adsorption [7][8][9][10][11][12], image reconstruction [13][14][15][16], tessellation [1,[17][18][19][20], and stochastic modeling [4,[21][22][23][24][25][26][27][28]. Among these methods, stochastic modeling has attracted attention for its utility in mimicking the real structure of FMs de novo. In stochastic modeling, all FM structural properties are drawn from various parameter distributions. Ideally, if those distributions are extracted from experimental measurements, the model should be empirically consistent. Despite the success of stochastic FM models, there remain challenges and issues that should be addressed to further their potential. Specifically, popularized stochastic methods for treating fiber orientation and tortuosity often lack tractability, transferability, or empirical consistency.
Many empirical studies analyzing the fiber orientation distributions of FMs have been reported. These characterization studies have conclusively shown that fibrous materials exhibit a range of material symmetries, such as nearly isotropic [29][30][31][32][33], transverse isotropic [26,[29][30][31][34][35][36][37][38][39], orthotropic [31][32][33]36,[39][40][41][42][43][44][45], and hexagonal [13,31,36,40,43,45]. Furthermore, fiber orientations often do not possess symmetry and are distributed anisotropically [32,[39][40][41]45]. Therefore, ideally, a stochastic FM model should have a means to recreate these material symmetries and to incorporate asymmetries. Although popular distributions such as those reported by Xia [46][47][48][49][50] and Schladitz [21][22][23]27,51] are undoubtedly useful for stochastic modeling and have been successfully used to mimic the orientation distributions of FMs, they lack transferability because they are specially formulated for specific symmetries. Fortunately, distributions that can model arbitrary material symmetries and asymmetries have been rigorously studied and employed in neuroscience to map fiber orientations within neurological tissue. Prominent distributions that possess a high level of transferability and empirical consistency include the Watson [52,53] and de la Vallee Poussin [52][53][54][55] distributions. These distributions possess three properties that allow them to be seamlessly employed to recreate arbitrary fiber orientation distributions: First, the distributions are antipodally symmetric; that is, the fiber orientationn f is equivalent to the fiber orientation −n f -a desirable simplification when modeling fiber orientation distributions. Second, multiple instances of the distributions can be superimposed without algebraic difficulty, which is the critical property giving rise to their flexibility. Lastly, these distributions have been thoroughly studied and translated to a vector form (convenient for computational work), a polar/azimuthal angle form (convenient for data analysis), and a tensor form. Moakher and Basser have published a comprehensive review on these distributions, where all of these forms are thoroughly discussed and documented [52]. However, regardless of a distribution's transferability, it cannot be directly used to accurately model 3D FM structures using 2D data extracted experimentally without correcting for projection distortion. Projection distortion occurs when 3D orientation distributions are projected onto a plane to be reported as a histogram or a polar plot. This is the same phenomenon that distorts the globe when projected onto a map. Without a relationship to correct for the projection distortion, fitting 2D orientation data with any probability distribution function (e.g., Xia, Schladitz, Watson, de la Vallee Poussin) does not yield proper parameters for recreating the 3D FM structure. Lastly, mitigating fiber penetration is often a necessary task in stochastic FM modeling, and may significantly affect the final fiber orientations, but the extent to which fiber orientations are disturbed during this process has not been rigorously quantified.
In addition, the treatment of curvilinear fibers in stochastic FM modeling is currently incomplete. The complex shapes that a fiber traces through space are typically reduced to a scalar quantity such as tortuosity [56][57][58][59] or straightness [40,[60][61][62]. Methods based on random walks and perturbations have been useful for stochastically modeling tortuous fibers [21][22][23]25,63,64]. However, these popular methods lack a connection between the input parameters of the random walk/perturbation and the tortuosity of the generated path. Without this relationship, attempting to model fiber tortuosity drawn from any distribution is a matter of trial and error. For this reason, the current treatments of fiber tortuosity in stochastic modeling are not tractable and, hence, cannot be empirically consistent until this issue is addressed. Moreover, as is an issue with fiber orientation, fiber tortuosity is subject to change when non-penetration conditions are enforced in FM modeling. The extent to which the non-penetration conditions affect the fiber tortuosity remains to be quantified.
Here, pertinent issues in stochastic FM modeling are addressed. First, the issue of projection distortion is addressed by developing a novel relationship to map from 2D fiber orientation data to 3D de la Vallee Poussin orientation distributions. The corrections from this mapping are shown to reduce parameter errors by up to 100%. Second, a novel empirical relationship is developed between the von Mises-Fisher random walk and the tortuosity of the path that it generates. The relationships for selecting von Mises-Fisher random walk parameters are shown to speed up the simulation of tortuous fibers bỹ 1200-fold and eliminate the issue of estimating parameters. Third, the effects of enforcing non-penetration conditions on fiber length, orientation, and tortuosity distribution are reported. Furthermore, the pore size of the generated FMs resulting from this model is measured and validated against known empirical relationships. Consequently, multivariate predictors for the mean pore size as a function of the FM's structural parameters are revealed. The results are consistent with the known inverse power law for mean pore size versus volume fraction and the linear scaling of mean pore size with mean fiber diameter [20,63,[65][66][67][68].

Methods
As discussed in Section 1, fibrous material structures can be modelled using a handful of random processes to reproduce the distribution of fiber diameter, length, orientation, and tortuosity observed in real structures. Here, the methods for modelling fiber diameter, length, orientation, and tortuosity are detailed. The methods employed to enforce nonpenetration conditions and measure the pore size of the resulting models are also discussed.

Model Domain, Volume Fraction, and Spatial Distribution of Fibers
Modeling began by declaring a spatial domain Ω, the fiber volume fraction φ f , mean fiber length l f , and mean fiber diameter d f . The spatial domain was a box whose faces met at right angles with side lengths L Ω and volume V Ω = L 3 Ω . The number of fibers, N f , was estimated as N f ≈ 4V Ω φ f π l f d f 2 , and the result was rounded to the nearest whole number. Then, N f seeds, x 0 , were uniformly distributed in the domain x 0,i = L Ω [N ur,1 , N ur,2 , N ur,3 ], where N ur,i is a uniformly random number between 0 and 1.

Drawing Random Numbers from Probability Distributions
In FM stochastic modeling, the geometric parameters of individual fibers-such as fiber orientation and tortuosity-are randomly drawn from distributions. Random numbers were drawn from a continuous probability distribution p(X) or histogram data H(X) using rejection sampling, where X is a randomly distributed variable. First, a number was uniformly sampled between [a, b] using X = (b − a)N ur + a, where N ur is a uniformly random number between 0 and 1, X is a sample value, and a and b are the lower and upper bound of the sampling domain, respectively. A second number, p , was uniformly sampled between [0, max(p(X))] for probability densities or [0, max(H(X))] for histograms. If p < p(X ), the sample value X was accepted, and if p > p(X ), the sample value X was rejected. When using histogram data, H(X ) was calculated by linear interpolation. For continuous probability distributions with infinite support, a was taken as µ − 6σ and b was taken as µ + 6σ, where µ and σ are the mean and standard deviation of the distribution, respectively. For continuous probability distributions with lower or upper finite support, a was set equal to the lower support and b was set equal to the upper support. For histogram data, a was taken as the left edge of the leftmost bin, and b was taken as the right edge of the rightmost bin. For 3D directional distributions, the same algorithm applied except that the distributed variable was a unit vector sampled from the unit sphere.

Fiber Diameters and Lengths
In our fiber model, the fiber diameters had a mean d f and a standard deviation std d f , while the fiber lengths had a mean l f and a standard deviation std l f . Both fiber lengths and diameters were drawn from a gamma distribution: where p Γ (X) is the probability density for a fiber with diameter or length X, k Γ is the shape parameter, θ Γ is the scale parameter, and Γ(k Γ ) is the gamma function evaluated at k Γ . The gamma distribution has a lower finite support at X = 0 and an infinite upper support. The parameters for the gamma distribution were calculated from the mean and standard deviation of X using the following equations:

Fiber Orientations
If a fiber is treated as a straight line, then a vector tangent to that line is defined as the fiber orientation. For a curvilinear fiber, the average of all tangent vectors along its backbone defines its mean orientation. A distribution of fiber orientations shows the material's symmetry, such as transverse isotropic (Figure 1(a.i)), orthotropic (Figure 1(a.ii)), or hexagonal (Figure 1(a.iii)). Here, fiber orientations were drawn from the de la Vallee Poussin (VP) probability distribution: where p VP n f is the probability density for a fiber (Figure 1(a.i-a.iv)) with orientationn f , unit vectorm VP is the preferred or mean orientation, and k 3D VP is the VP orientation strength parameter. Equation (4) is applicable to describe transverse isotropic fiber orientation distributions when k 3D VP > 0 or isotropic fiber orientation distributions when k 3D VP = 0. Any arbitrary orientation distribution can be treated as a linear combination of multiple VP distributions of different kernels: where N VP is the number of superimposed VP kernels. For isotropic and transverse isotropic orientation distributions N VP = 1; orthotropic N VP = 2; hexagonal N VP = 3; and asymmetric N VP can be arbitrary. Equations (4) and (5) were fitted to 2D experimental histogram data using the polar (θ) and azimuthal (φ) angle forms of the VP distribution. To do so, letn f = (cosθsinφ, sinθsinφ, cosφ) andm VP,i = (cosθ VP,i sinφ VP,i , sinθ VP,i sinφ VP,i , cosφ VP,i ) with 0 < θ, θ VP < 2π and 0 < φ, φ VP < 2π. The mathematics can be simplified by setting φ = φ VP = π/2, and then Equations (4) and (5) can be rewritten as follows: Polymers 2022, 14, x FOR PEER REVIEW 5 of 21

Fiber Tortuosity
The idealized picture of a fiber is a straight cylinder; a slightly more realistic description is a curvilinear cylinder. In practice, there are many measurements to define the curviness of a fiber or a path, such as straightness, sinuosity, fractal dimension, or tortuosity, with tortuosity being one of the most popular. Tortuosity, τ f , is the ratio of a curve's length to the total distance traveled and is always greater than or equal to one. In this paper, fiber tortuosity τ f was drawn from probability distributions or interpolated histograms and followed a lognormal distribution: where p ln τ f is the probability density of a fiber having a tortuosity τ f , while m ln and s ln are the expected value and standard deviation of τ f 's natural logarithm, respectively. Both m ln and s ln can be calculated from the mean and standard deviation of τ f , as follows: Tortuosity can never be less than 1, so Equation (8) is truncated at τ f = 1, meaning that Equations (9) and (10) are only approximations here.
Curves with tortuosity τ f ,i were generated using the von Mises-Fisher (VMF) directional probability distribution: to guide random walks that were subsequently smoothed with cubic splines (Figure 1(b.i,bii)).
Here, p V MF (n V MF ) is the probability of stepping in directionn V MF ,n f is the fiber orientation, and κ V MF is the VMF strength parameter. A fiber was generated from a smoothed VMF random walk with step size ∆ V MF , spline step size ∆ CS , and preferred directionm V MF =n f (Figure 1(b.i)). For a low orientation strength κ V MF , the simulated fibers were highly tortuous ( Figure 1(b.ii)), while for high orientation strength κ V MF they were less tortuous ( Figure 1(b.ii)). Random walks began at fiber seeds x 0 ; a direction was drawn from the VMF distribution (Equation (11)), and a step was taken using The number of steps taken was calculated as Afterwards, each random walk was smoothed with cubic spline interpolation with step size ∆ CS . The criteria for choosing κ V MF , N V MF , ∆ V MF , and ∆ CS are detailed in Section 3.2.

A Fiber as a Chain of Particles and Non-Penetration Conditions
After employing the methods discussed in Sections 2.1-2.5, the result was a box full of randomly oriented, randomly positioned fibers. In this case, it is practically inevitable that some or many fibers are penetrating one another. This penetration needs to be addressed in many cases, especially when simulating non-rigid FM models. Here, classical mechanics were used to mitigate fiber-fiber penetration by introducing a repulsive force between fibers. In attempt to maintain the prescribed fiber tortuosity and length, restoring forces for bond lengths and angles were included between particles in a fiber. Lastly, the non-penetrating configuration is always a lower energy state than the penetrating configuration. A drag force was introduced to equilibrate the system towards the lower-energy non-penetrating configuration. A single fiber was represented as a chain of spherical particles x f , whose positions were defined by the smoothed VMF walk. Periodic boundary conditions were used to reinsert particles initially existing outside the spatial domain. Harmonic bonds, cosine angles, and soft contact potentials were used to enforce non-penetration conditions (Figure 1(c.i-c.iii)), as follows: where f b , f θ , and f c are the forces felt by particles due to the bond, angle, and contact potential, respectively; k b , k θ , and k c are the stiffnesses of the bond, angle, and contact potential, respectively; r ij is the distance vector between particles i and j; R ij is the distance between particles i and j; R 0,ij is the equilibrium distance between particles i and j; R c,ij is the cutoff distance of the contact potential; θ ijk = acos r ij ·r kj /R ij R kj is the angle between particles i, j, and k; and θ 0,ijk is the equilibrium angle between particles i, j, and k. The equilibrium distance between two bonded particles in the fiber was taken as ∆ CS of the smoothed VMF walk. The equilibrium angle between three consecutive particles was taken as the initial angle between them after the smoothed VMF walk. The contact potential cutoff distance between particles i and j was taken as . The contact potential was not calculated between bonded particles or particles participating in a bond trio.
Non-penetration conditions were enforced by running damped classical mechanics simulations. A drag force was added to each particle using Stoke's law (Figure 1(c.iv)).
where η is the viscosity of a surrounding fluid and v i is the velocity of the particle. Using the forces defined in Equations (12)-(16), the velocities v and positions x of the particles were updated by the velocity Verlet time-stepping algorithm. The time step ∆t was taken as . Periodic boundary conditions were employed during the dynamics simulation, and fiber ends were held fixed for the first 10,000 timesteps to restrict fiber rotations during high-energy dynamics.
The fiber elastic modulus E, Poisson's ratio ν, particle mass m, and fluid viscosity η were chosen randomly, and linear elasticity theory was used to relate E to {k b , k θ , k c }. Bond stiffness, k b , was set equal to EA/R ij,0 , where A = π 4 d 2 f is the cross-sectional area of the fiber. Angle stiffness, k θ , was set equal to EI/R ij,0 , where I = π 32 d 4 f is the second moment of area of the fiber. Contact stiffness, k c , was set equal to R 9/4 ν 1/3 E/16, where 1/R = 1/d f ,i + 1/d f ,j (approximately close to Hertz contact). Hertzian and other contact theories were applicable; however, the soft potential in Equation (15) was chosen to avoid large accelerations when the initial separation of two particles was much smaller than the sum of their radii.
After constructing the FM networks (Figure 1d), their corresponding structure parameters (such as fiber volume fraction, orientation, length, tortuosity, etc.) were characterized and compared to the distributions initially used for each simulation.

Pore Size Measurement
Once a simulated network was finished, a watershed method [19,20,69] was used to calculate the pore diameters, d p (Figure 1e). First, the simulation domain was segmented into voxels (Figure 1e, top left). The voxel edge lengths were 5 times smaller than the smallest fiber diameter. If a voxel was contained or partially contained in a fiber particle, the voxel was given a value of 1; otherwise, it was given a value of 0 (Figure 1e, top right). The distance transform was applied to the voxels, where each voxel was assigned a value corresponding to the minimum distance from itself to a voxel with a value of 1 (Figure 1e, bottom left). By convention, the distance transform was negated, and the voxels located in local minima were marked as pores (Figure 1e, bottom right). Local minima were found using the watershed method. The pore diameters were taken as the value of the distance transform at local minima.

Relating de la Vallee Poussin Parameters Measured from 2D data to 3D Structures
Fiber orientation is the mean direction that a fiber points towards, calculated by averaging the tangent vectors along the axial direction of a fiber. The distribution of fiber orientations is typically modelled with a 3D distribution but measured and reported as a 2D distribution. In Figure 2a, two transverse isotropic VP distributions with the same orientation strength are shown (k 3D VP = 6). The blue VP distribution lies in the projection plane, while the orange VP distribution has been rotated out of the "projection plane" by an angle ψ = 45 • . An orientation that could have been drawn from either distribution is shown asn f . Since experimentally the observer is measuring the orientation from above (i.e., projecting the 3D network into a 2D network), the effective orientationn f being measured is the projection ofn f on the projection plane. Such a projection process could result in projection distortion, i.e., the observed 2D orientation distribution would vary when the ψ is changed with respect to the observer even though the 3D distribution shape remains the same. For example, two fiber orientation distributions with the same 3D fiber orientation strength (k 3D VP = 6) were generated according to the blue and orange distributions in Figure 2a, and the corresponding 2D orientation distributions on the horizontal projection plane were obtained. Figure 2b shows the corresponding orientation distribution histograms, and clearly the distribution for ψ = 0 • (blue histogram) is significantly different from that for ψ = 45 • (red histogram). Equation (6) was used to fit these two histograms, and the orientation strength parameter k 2D VP was extracted to be k 2D VP = 15 for ψ = 0 • and k 2D VP = 8 for ψ = 45 • . Note that the true 3D orientation strength k 3D VP = 6 is the same for both cases. Hence, without understanding how the projected 2D measurement affects the fiber orientation distribution, the "true" distribution cannot be recovered from 2D data.
In order to obtain the relationship of the 2D and 3D orientation distribution, a detailed investigation was conducted by systematically varying both k 3D VP (from 0 to 10) and ψ (between 0 • and 90 • ). Figure 2c plots the extracted k 2D VP as a function of k 3D VP for different ψ for transverse isotropic distribution. At a fixed ψ, the measured k 2D VP increases linearly with k 3D VP , while the slope of the k 2D VP vs. k 3D VP plot changes as a function of cos 2 ψ (Figure 2d; filled circles are data obtained from Figure 2c, and the solid curve is the fitting). Similar measurements and analyses were conducted for orthotropic and hexagonal distributions and showed identical trends. In fact, our systematic investigation revealed the following relationship: According to Equation (18), if one knows the relative orientation between the 3D orientation distribution and the projection plane of the measurement, one should be able to extract the true 3D orientation strength based on the measured value. Alternatively, if experimentally multiple orientation samplings were obtained, then: represents the experimental fiber orientation data; "VP fit" is the regression of Equation (7) to experimental data; the orange histogram is the fiber orientation distribution of simulated FM structures when projection distortion is accounted for; the white histogram is the fiber orientation distribution of simulated FM structures when directly using k 2D VP (from regression) as the orientation strength parameter k 3D VP (for drawing random fiber orientations).
According to Equation (18), if one knows the relative orientation between the 3D orientation distribution and the projection plane of the measurement, one should be able to extract the true 3D orientation strength k 3D VP based on the measured k 2D VP value. Alternatively, if experimentally multiple orientation samplings were obtained, then: Thus, the true k 3D VP can be obtained experimentally by measuring the average k 2D VP at different ψ orientation sampling. If the limits of integration span an integer multiple of π (Ψ 2 − Ψ 1 = nπ), then Equation (19) can be reduced to k 2D . To further illustrate the relationship (Equation (18)) between the 3D and 2D orientation distributions, we simulated the FMs fibrin [29], stone wool [13], cross-laid polyester [43], and collagen [41] using experimental orientation data. The discrete data shown in Figure 2(e.i-h.i) are the experimentally measured 2D fiber orientation distributions. These data were fitted by the 2D VP distribution (Equation (7)); see the solid blue curve in each figure), and the corresponding fitting parameters are summarized in Table 1. For fibrin, stone wool, cross-laid polyester, and collagen, we obtained k 2D VP = 0.98, k 2D VP,1 = k 2D VP,2 = k 2D VP,3 = 5.25, k 2D VP,1 = k 2D VP,2 = 2.5, and [k 2D VP,1 , k 2D VP,2 ] = [4, 2.2], respectively. Since experimentally no ψ-dependent measurements were conducted, we assumed that that ψ = 0 • for all four cases. Thus, based on Equation (18), the corresponding 3D orientation strength k 3D VP can be estimated as shown in Table 1. Visually it is quite difficult to tell the difference between the four FMs. The projected 2D orientation distributions were obtained from the simulated FMs and are plotted as orange histograms in Figure 2(e.i-h.i). Visually, the 2D orientation distributions obtained from the simulated FMs match quite well with the experimental distributions, with a slightly larger variation for collagen (Figure 2(h.i)). Moreover, as a comparison, when we assumed that k 3D VP = k 2D VP (which has been often done in experiments, and we call it non-corrected situation), another set of FM models were generated for the above four cases, and the corresponding 2D orientations are plotted as the white histograms in Figure 2(e.i-h.i). Compared to the corrected projection distortion situation-i.e., the orange histograms-the non-corrected situation generally overestimates k 3D VP . Therefore, the experimentally measured 2D orientation distribution parameters cannot be directly taken as parameters for 3D orientation distributions as usual, because this can lead to errors of up to 100% when estimating k 3D VP . A good approximation is to use Equation (18) to convert the 2D parameters to 3D parameters. Based on both the experimental and fitting results for the 2D orientation distributions, the orientation distributions of fibrin, stone wool, and cross-laid polyester are quite symmetric, while that of collagen is not symmetric. Clearly, even with the distortion correction, the obtained [k 2D VP,1 , k 2D VP,2 ] from the simulations does not match well with those from the experimental data of the asymmetric case.

Relating κ V MF and N V MF to Fiber Tortuosity τ f
Fiber tortuosity is the ratio of a fiber's contour length to its end-to-end length (Figure 3a). Controlling fiber tortuosity can be a tricky task in simulations because it is difficult to know the end-to-end distance of randomly generated curves in advance. This section details the analysis of random walks guided by the VMF probability density and smoothed with a cubic spline. distributions, the orientation distributions of fibrin, stone wool, and cross-laid polyester are quite symmetric, while that of collagen is not symmetric. Clearly, even with the distortion correction, the obtained [ , , , ] from the simulations does not match well with those from the experimental data of the asymmetric case.

Relating and to Fiber Tortuosity
Fiber tortuosity is the ratio of a fiber's contour length to its end-to-end length ( Figure  3a). Controlling fiber tortuosity can be a tricky task in simulations because it is difficult to know the end-to-end distance of randomly generated curves in advance. This section details the analysis of random walks guided by the VMF probability density and smoothed with a cubic spline. represents the experimental data for the FMs collagen [70], sintered metal fiber [56], and fiberboard [58]; τ is the measured fiber tortuosity distribution from the simulated FMs. represents the experimental data for the FMs collagen [70], sintered metal fiber [56], and fiberboard [58]; τ sim f is the measured fiber tortuosity distribution from the simulated FMs.
The VMF random walk has four independent variables: the orientation strength parameter κ V MF , the number of random steps N V MF , the step size ∆ V MF , and the separation of particles in the cubic spline ∆ CS . The length of the smoothed random walk is approximately equal to N V MF ∆ VFM . Random walks with orientation strengths between 0.001 and 300, a number of steps between 3 and 80, step sizes between 3 and l f 800 were used to generate 5800 fibers to study the relationships between smoothed random walk parameters and fiber tortuosity. Figure 3a shows three representative fibers with different τ f ; when τ f increases from 1.05 to 1.35, the fiber changes from a fairly straight fiber to a curved fiber with multiple bends. Preliminary analysis showed that ∆ CS had little-to-no effect on the tortuosity of a path generated by the VMF random walk. However, κ V MF strongly influenced the tortuosity of the resulting path, and N V MF limited the tortuosity that a random walk could achieve. The variables N V MF and ∆ V MF were correlated with one another; thus, only one of them needed to be studied here. The number of VMF steps was chosen as a control variable for the response variable τ f . In Figure 3b, the control variable κ V MF is plotted on the y-axis (log), and the dependent variable-the expectation value of tortuosity τ f -is plotted on the x-axis, as is customary when inverting a problem. As is already known, the expected fiber tortuosity τ f was inversely proportional to the orientation strength parameter κ V MF . However, decreasing κ V MF did not necessarily increase τ f . If the number of VMF steps N V MF was too small, decreasing κ V MF led to no change in τ f . This can be observed in Figure 3b,c (solid black curves); as κ V MF decreases, τ f approaches a limit value, τ f limit , for a fixed value of N V MF . The orientation strength parameter was found as a function of the expected fiber tortuosity: given that the tortuosity of the walk was not limited by the number of steps (black curve in Figure 3b). The observed relationship was the exponential of a rational function with a vertical and slant asymptote (Equation (20)). The observation that the expected fiber tortuosity can be limited by the number of steps in the random walk required a second equation be defined for determining an appropriate number of steps, N V MF . In Figure 3c, the minimum number of steps required for producing curved fibers with the expected tortuosity τ f is shown (black curve). The region below the black curve in Figure 3c is inaccessible, i.e., the random walk with N V MF steps that produces curves with τ f does not exist. The minimum number of VMF steps, N * V MF , required to produce curves with expectation value τ f is given as follows: The colored dashed curves in Figure 3c represent N V MF using different multiples of N * V MF via: where β is the number multiplied by Equation (20) to determine the number of steps taken (Figure 3c). Using more steps than necessary, β > 1, to produce curves with mean tortuosity of τ f was found to reduce the standard deviation of fiber tortuosity for a fixed κ V MF . For instance, in Figure 3c, for β = 5 the mean absolute percentage error (MAPE) between the tortuosity τ f for 1000 simulated fibers and the expected tortuosity τ f is ∼ 20%. When β is increased to 20, the MAPE decreases to ∼ 20%. Therefore, as N V MF approaches infinity, a single value of κ vm f maps to a smaller range of possible fiber tortuosity τ f . If the goal of the random walk is to produce a curve whose tortuosity is very close to τ f , then increasing β will reduce the mean absolute percentage error (MAPE) of that procedure (Figure 3c). Since, for a finite number of steps, a single value of κ V MF maps to a range of possible τ f , the curves can be screened for error. A demonstration of the simulated FM models using smoothed VMF random walks and criteria given by Equations (20)- (22) with β = 20 is shown in Figure 3(d.i-f.ii), based on the experimental fiber tortuosity distributions for collagen [70], sintered metal fiber [56], and fiberboard [58]. The scattered data in Figure 3(d.i-f.i) are the experimental fiber tortuosity distributions, which are directly used to randomly select the fiber tortuosities for simulated networks (Figure 3(d.ii-f.ii)). Unlike other fiber properties, such as fiber length or diameter, a general probability density function for fiber tortuosity is elusive. Therefore, it is preferable to directly select fiber tortuosity from experimental histogram data when they are available, rather than from fitted probability densities. Since Equation (20) only defines κ V MF in terms of the expected fiber tortuosity, the process is not deterministic. To increase accuracy further, the random walks were screened. If the error in the tortuosity of a generated curve was greater than 5% of the prescribed tortuosity, the walk was discarded. Due to cubic spline smoothing, the lengths of fibers generated by the random walks were not exactly equal to N V MF ∆ V MF . Some representative FM networks are shown in Figure 3(d.ii-f.ii). Since the τ f for collagen is the smallest, the fibers shown in Figure 3(d.ii) are quite straight, whereas the τ f for fiberboard is maximized, and all of the fibers in Figure 3(f.ii) are very curvy. The extracted tortuosity distributions from the simulated FM networks are plotted as orange histograms in Figure 3(d.i-f.i), and they match quite well with the experimental data. The MAPEs between the randomly selected fiber lengths and tortuosities and the realized fiber lengths and tortuosities are reported in each plot, and they are quite small. Such a result demonstrates that using the strategy proposed in Section 2.5 we can successfully incorporate tortuosity and tortuosity distribution in FM simulation.
A set of simulations were conducted to compare the selection of VMF random walk parameters from Equations (20)- (22) with selecting them randomly. We generated 1000 fibers to reproduce the sintered metal fiber tortuosity distribution (tortuosity between 1 and 3.5). A range of 10 ≤ N V MF ≤ 100 and 0.1 ≤ κ V MF ≤ 100 was used to recreate the situation of a naïve guess and check. For comparison, Equations (20)- (22) with β = 20 were used to generate a statistically identical set of fibers. Fibers were only accepted in both cases if the absolute percentage error was less than 5%. For the naïve guess and check, the process took 40 min to complete, and the rejection rate was 986 fibers per 1000 randomly generated fibers. With Equations (20)- (22) and β = 20, the process was completed in 2 s, with a rejection rate of 32 fibers per 100 randomly generated fibers. This corresponds to a 1200-fold increase in speed, and Equations (20)-(22) remove the worry of guessing an appropriate range for κ V MF and N V MF .

Effect of Enforcing Non-Penetration on τ f ,n f , and l f Distributions
When objects with finite size are randomly placed in a domain, such as in stochastic FM modeling, there is some probability that they overlap or penetrate one another (Figure 4(a.i,aii)). This probability increases as the volume fraction increases. For medium and long fibers, volume fractions of~1% are sufficient to ensure that fiber overlap occurs. Such an effect may be undesirable in FM simulations. A typical method for treating fiber overlap is to treat the fibers as elastic objects that repel one another, as elaborated in Section 2.6. Polymers 2022, 14, x FOR PEER REVIEW 14 of 21    Figure 4(a.i)). After enforcing non-penetration conditions, the fiber overlap reduced substantially, and the fibers translated, rotated, and deformed (Figure 4(a.ii)). Figure 4b shows the change in percent fiber overlap (PFO) during the simulation. Percent fiber overlap was calculated as the total overlapping volume divided by the total fiber volume multiplied by 100. It is clear that enforcing non-penetration conditions by treating fibers as elastic bodies that repel one another is an effective means to reduce fiber penetration. However, it is also the case that this procedure can change fibers' shapes and orientations. Therefore, the tedious work of matching fiber length, tortuosity, and orientation distributions is threatened by this procedure. Although no thorough remedy can be offered here, we statistically quantified the effects of enforcing non-penetration conditions on the fiber shape and orientation changes as well as the overall FM model. Figure 4c shows the distribution of changes in fiber lengths, tortuosities, and orientations as a result of enforcing non-penetration conditions. The boxplots were generated from 100 fiber network simulations with mean fiber aspect ratios l f /d f between 10 and 150, fiber volume fractions φ f between 1% and 40%, fiber orientation strengths k 3D VP between 0 and 10, and mean fiber tortuosities τ f between 1.05 and 2.5, and included FM models with isotropic, transverse isotropic, orthotropic, and hexagonal symmetries. Individual fiber lengths, tortuosities, and orientations changed by~2%,~14%, and~5% on average. However, the changes to the mean fiber length, tortuosity, and orientation distributions were smaller (Figure 4c):~1%,~9%, and~3%, respectively. While some fibers became longer, others became shorter; while some fibers became straighter, others became curvier; and while some fibers rotated clockwise, others rotated counterclockwise. The data provided here show that enforcing non-penetration conditions has a significant effect on the final structure of the FM models-most notably on fiber tortuosity and orientation. These results suggest that more work is needed to develop better criteria and methods for enforcing non-penetration conditions-particularly criteria and methods aimed at preserving fiber tortuosity and orientation. Until then, these data serve as a rough estimate for the unavoidable errors generated in stochastic FM modeling when enforcing non-penetration conditions. This section is intended to provide actions that may be taken to reduce the errors (for fiber length, orientation, and tortuosity) that occur when employing non-penetration conditions. However, the complexities of this issue are extreme enough to warrant a dedicated study of the topic, and no quantitative relationships can be confidently offered here. That being said, the magnitude of errors generated by enforcing non-penetration conditions is significant, as briefly discussed above.

Pore size Distribution Analysis of FM Models
FMs are inherently porous, and the pore structure is a direct consequence of the FM structure. Pore size d p is known to increase linearly with increasing fiber diameter d f and to decrease with the square root of the volume fraction φ f [66,67,[71][72][73][74][75]. A final step in verifying the structure of an FM model is to quantify the pore structure and compare it to empirical trends. The mean pore size d p of 850 FM models was measured. The models had volume fractions φ f between 0.01 and 0.65, fiber diameters d f from 1.5 nm to 4 mm, fiber lengths l f from 25 nm to 400 mm, fiber tortuosities τ f from 1.03 to 6, and fiber orientation strengths k 3D VP from 0 to 20. Univariate (Figure 5a-c) and multivariate (Figure 5d-f) analyses were conducted using traditional linear regression to relate d p or pore size distributions to d f , l f , φ f , τ f , and k 3D VP . Figure 5a shows the empirically consistent result with a linear relationship d p ∝ d f [65][66][67][68][71][72][73][74], while Figure 5b shows the empirically consistent result where [65][66][67][68][71][72][73][74]. The data show little-to-no relationship of d p with l f , τ f , and k 3D VP . As shown in Figure 5c, d p appears to decrease linearly with τ f . The effect is small compared to the relationships of d p vs. d f and d p vs. φ f , especially over typical ranges of fiber tortuosity (~1-2).
Knowing the mean fiber diameter, the mean pore size could be estimated within ~60% error ( Figure 5(d.ii)) using Equation (23). Volume fraction is the second most significant predictor, as shown in Figure 5(c.i), with a semi-log plot. The predicted relationship is as follows: Knowing the mean fiber diameter and volume fraction of the FM allows for estimations of mean pore diameter within ~40% error ( Figure 5(e.i,e.ii)) using Equation (24). Lastly, knowing the mean fiber diameter, tortuosity, and volume fraction, the mean pore diameter could be estimated within ~20% error ( Figure 5(fi-e.ii)) using Equation (25): Multivariate analysis would be more insightful. Figure 5(d.i) shows that the predominant predictive variable of mean pore diameter is mean fiber diameter, with the following relationship: Knowing the mean fiber diameter, the mean pore size could be estimated within~60% error ( Figure 5(d.ii)) using Equation (23).
Volume fraction is the second most significant predictor, as shown in Figure 5(c.i), with a semi-log plot. The predicted relationship is as follows: Knowing the mean fiber diameter and volume fraction of the FM allows for estimations of mean pore diameter within~40% error ( Figure 5(e.i,e.ii)) using Equation (24). Lastly, knowing the mean fiber diameter, tortuosity, and volume fraction, the mean pore diameter could be estimated within~20% error ( Figure 5(fi-e.ii)) using Equation (25): Reasonably speaking, since very high fiber tortuosities were included in the dataset, the error in predicting the mean pore diameter when knowing only the mean fiber diameter and volume fraction would likely be less than what reported here in practice. On the other hand, Equation (23) may be inaccurate when predicting mean pore diameter for very lowor high-volume-fraction FMs. Equations (23) and (24) report well-known pore diameter relationships and support the validity of stochastic FM modeling. However, Equation (25) suggests that experimental investigation of the effect of fiber tortuosity on pore structure might be warranted.

Limitations
In Section 3.1, a relationship for correcting de la Vallee Poussin orientation probability density parameters is presented. It shows that the prescribed corrections are useful for transverse isotropic, orthotropic, and hexagonal fiber orientation symmetries. However, the prescribed corrections are not useful for asymmetric fiber orientations. Furthermore, unless multiple measurements are taken, one must assume that due diligence has been taken to approximately align the FM's primary fiber orientation directions with the imaging plane. In Section 3.2, relationships are provided to select parameters for a VMF random walk to produce curves with the desired tortuosity. The relationships perform well for their intended purpose; however, random walks are not deterministic; therefore, the tortuosity of a random walk cannot be predicted deterministically. The provided relationships help in selecting parameters that correspond to a high probability of producing a random walk with a desired tortuosity. Therefore, screening of random walks is still advised. In Section 3.3, the magnitude of errors that arise from enforcing non-penetration conditions is quantified. This section shows the limitations of stochastic FM modeling in cases where fiber overlap is not permissible in an FM model.

Conclusions
Stochastic modeling is a useful tool for studying FMs with numerical methods. This report addresses three critical issues in stochastic FM modeling methods, namely, projection distortion in fiber orientations, fiber tortuosity, and non-penetration conditions. Analysis of fiber orientation distributions revealed that 2D histograms of fiber orientation distributions suffer from projection distortion. This results in mismatches between the fitted orientation distributions and reality. We provided a relationship between the 2D and 3D de la Vallee Poussin probability distribution functions to correct for projection distortion. The relationship was valid and accurate for transverse isotropic, orthotropic, and hexagonal material symmetries, but inaccurate for anisotropic distributions. Further work is needed to extend the relationships to arbitrary material symmetries. After analyzing the tortuosity of paths generated by VMF random walks, we were able to develop criteria to relate the VMF random walk parameters to the expected tortuosity of the path. It is known that enforcing non-penetration conditions results in altering the shape and orientation of fibers in stochastic FM models. The extent to which fiber lengths, tortuosities, and orientations change has not been studied and reported in the literature. Here, we showed that fiber lengths changed by~2% on average, fiber orientations changed by~5% on average, and fiber tortuosities changed by~14% on average, while the mean fiber length changed bỹ 1% on average, the mean fiber orientation strength changed by~3% on average, and the mean fiber tortuosity changed by~9% on average. While the change in fiber lengths was quite low, the change in fiber tortuosities and orientations was very significant for the cases studied here. Further work is needed to develop methods and criteria for enforcing non-penetration conditions, which better preserve fiber tortuosities and orientations. Lastly, an analysis of the pore structures produced by these procedures was conducted. The analysis showed that the mean pore size of these models scaled linearly with fiber diameter and was inversely proportional to the square root of the volume fraction. These results are consistent with known empirical trends. Furthermore, the analysis showed that fiber orientations and lengths had little-to-no effect on the mean pore diameter, but fiber tortuosity showed a meaningful effect. Specifically, large fiber tortuosities (greater than 2) significantly changed the mean pore diameter. Since fibers in real FMs are typically not that large, this effect is functionally negligible. However, this result may warrant experimental investigation. Further work is necessary to develop a relationship for correcting orientation distribution parameters in cases where fiber orientations are asymmetric and for predicting errors that occur when enforcing non-penetration conditions so that they can be dealt with preemptively.
Author Contributions: Y.Z. and X.W. conceived the idea. N.F. developed the theory, performed the computations, and analyzed the data. All authors have read and agreed to the published version of the manuscript.