Semi-flexible polymer networks underlie the mechanical behavior of a wide range of natural and engineered materials. Interestingly, these networks are often embedded within amorphous matrices rather than appearing in isolation. However, despite their frequent occurrence as embedded rather than isolated networks, few prior studies have specifically investigated the role of embedding on the emergent mechanical behavior of these systems. To address this, we adopt a mortar-type embedding approach within the finite element framework and perform simulations to systematically fill this knowledge gap. Within this study, we focus on soft tissues as an exemplary class of materials where embedded semi-flexible polymer networks are essential to mechanical function. Specifically, we investigate the role of embedding on the strain energy distribution within the networks across the bending, stretching, torsional, and shear fiber-level loading modes. In addition to revealing the role of embedding on the networks themselves, we also investigate how the networks affect the mechanics of the matrix material. Together, we find that embedding fundamentally alters the mechanics of the networks and the surrounding matrices. Most importantly, we find that embedding semi-flexible polymer networks leads to strain-stiffening and negative Poynting effect of the resulting composite material. Furthermore, semi-flexible polymer networks induce stress heterogeneity in their host material and increase its resistance to compression. Overall, our work improves our fundamental understanding of an important class of materials. By making our implementation openly available, we also hope to help others learn more about embedded semi-flexible polymer networks in the context of materials other than soft tissues. 1 Introduction Semi-flexible polymer networks are an interesting and important class of materials. They play critical roles in the textile industry, as filter materials, in energy storage devices, and many other engineering applications [ 1,2]. However, their most important role may be as the ubiquitous mechanical building blocks of macroscale organisms [ 3,4,5]. For example, semi-flexible collagen networks form the structural backbone of our bodies’ soft tissues, such as ligaments, skin, heart valves, blood vessels, and many more [ 6,7,8]. Similarly, semi-flexible fibrin networks form the structural backbone of blood clots [ 9,10], while actin networks give structure and (contractile) function to our cells [ 11]. Clearly, characterizing and modeling the mechanics of semi-flexible polymer networks is critical to our understanding of many engineering applications, but especially to our understanding of human health and disease. However, understanding the structure-function relationships of these building blocks is not without challenges. Critically, semi-flexible polymer networks exhibit – and are said to endow their host matrices, such as cells and tissues – with complex and emergent mechanical behaviors including strain stiffening and negative Poynting effect [12, 13, 14]. Mechanics of Embedded Fiber Networks We are not the first to recognize the need for characterizing and modeling the mechanics of semi-flexible polymer networks [ 15,16]. Of the many previous investigations, most have taken one of two approaches. They either homogenized the semi-flexible polymer networks using a generalized structural tensor or the angular integration method [17,18,19,20,21], or they modeled semi-flexible polymer networks discretely. That is, fibers and their interactions were discretely accounted for as 1D elements. For example, Leng et al, used 1D hyperelastic trusses to approximate the mechanics of polymer fibers [ 22]; others used Timoshenko beams [ 13,4,14,23,24,25], while others modeled the fibers as either tension-only elements or as axial springs [26, 27, 28]. Broadly speaking, we are motivated to contribute to this area of research because the vast majority of efforts to model semi-flexible polymer networks to date have considered these networks in isolation. However, both in many engineering applications and in nature, networks often do not appear in isolation. Instead, these networks are often embedded within amorphous matrices. For example, in soft tissues, collagen is embedded in what is often called the “amorphous ground substance.” [ 29] Similarly, in blood clots, fibrin is embedded within a sea of red blood cells, which may also be viewed as an amorphous embedding medium [ 30]. Relatively few prior efforts have attempted to directly model the interactions that arise from embedded semi-flexible biopolymer networks. Among them, Zhang et al. embedded discrete 1D fibers in a mesh of tetrahedral solid elements [ 31]. They discretized the domain so that the solid element edges coincided with the center lines of the fiber elements. Thus, the solid and the beam elements shared the nodal degrees of freedom. The same model was used in a later study of the effect of fiber crimp on the stress distribution of the solid matrix [ 32]. While important contributions, this and other similar approaches suffered from significant shortcomings. Specifically, they make the solution necessarily mesh-dependent. Additionally, these approaches limit the choice of fiber shapes and can be quite computationally expensive as they require a high mesh density around the discrete fibers to resolve the interface. Given the small number of prior efforts to characterize and model the mechanics of embedded semi-flexible polymer networks and their limitations, there are significant knowledge gaps. Critically, it is unclear if embedding is critical to understanding the fundamental mechanics of these systems. Thus, our goal in this paper is to fill this knowledge gap. To this end, we model fibers as connected Timoshenko beams and embed them within a hyperelastic matrix. We use the mortar-type finite element approach by Steinbrecher et al. that overcomes the need for conforming meshes at their interface [ 33]. Thereby, our models become computationally tractable while simultaneously allowing for exploration of complex fiber shapes and fiber-network geometries. The manuscript is organized as follows. In Section 2, we first briefly recall and verify our implementation of Steinbrecher et al.’s embedding technique. In Section 3, we show a single fiber example where we get a first glimpse of the effect of an embedding matrix. In Section 4, we conduct a detailed numerical analysis of semi-flexible polymer networks within an embedding matrix. Therein, we focus on two mechanical phenomena that are often attributed to semi-flexible polymer networks: strain stiffening and the negative Poynting effect. Finally, in Section 5, we conclude our work with a summary and a future outlook. 2 Computational Embedding Approach 2.1 Brief Summary of Steinbrecher et al.’s Numerical Approach We model semi-flexible polymer networks – subsequently referred to as “fiber networks” – as connected, discrete, spatial Timoshenko beams that are embedded in an isotropic, hyperelastic solid material, i.e., “matrix”. We then enforce the fiber-to-matrix coupling using a mortar-type finite element method, originally developed by Steinbrecher et al. [ 33]. This technique essentially constrains the fiber deformation field ufto that of the matrix um. It does so through a Lagrange multiplier field λthat may be interpreted as a line traction along the fiber opposing the relative displacement between fiber and matrix. Below we briefly describe essential components of the approach. We begin by recalling the spatial discretization of the matrix displacement field, the fiber displacement field, and the Lagrange multiplier field as 2 Mechanics of Embedded Fiber Networks um h=nmX k=1Nk(ξm, ζm, ηm)dm k (1) uf h=nfX l=1Φl(ξf)df l(2) λh=nλX j=1Φj(ξf)λj (3) where Nkare linear shape functions for the matrix nodes k,Φl/jare standard Lagrange shape functions (either linear or quadratic), and dm kanddf land are the nodal displacements of the matrix and fiber, respectively, while λjare the nodal values of the Lagrange multiplier field. These approximations give rise to the traditional stiffness matrices of the matrix domain Kmand the fiber domain Kf, as well as three additional (coupling) stiffness matrices as defined over the fiber domain Γf h: D(j,l)=Z Γf hΦjΦldsI3×3∈R3×3(4) M(j,k)=Z Γf hΦjNkdsI3×3∈R3×3(5) κ(j,j)=Z Γf hΦjdsI3×3∈R3×3(6) In turn, these coupling matrices define the relative displacement field and the Lagrange multiplier field as gc(dm,df) = [−M D ]dm df ∈Rnf×3, λ=ϵκ−1gc(dm,df)∈RnB×3, (7) where ϵis the penalty parameter that enforces the Lagrange multiplier field and thereby minimizes the relative displacement between matrix and fiber, i.e., enforces the embedding [ 33]. Note, integrating the matrices M,D, andκ poses several practical challenges. For example, integration domains may be overlapping as some beams will exceed the boundaries of one or even multiple matrix elements. Thus, selecting those degrees of freedom between the matrix and the fiber domain that are coupled is non-trivial and must be considered carefully during integration. Please see Supplement S1 and Fig. S1 for additional implementation details. Finally, we use a Newton-type iterative strategy to solve for the matrix and fiber degrees of freedom under the constraint of their coupling by solving the linearized problem:  Km+ϵMTκ−1M −ϵMTκ−1D −ϵDTκ−1M Kf+ϵDTκ−1D∆dm ∆df = fm+ϵMTκ−1gc ff−ϵDTκ−1gc , (8) where fmandffare the global force vectors for the matrix and the fiber nodes, respectively. We implement this approach through a user element (UEL subroutine) in Abaqus Standard (Abaqus/CAE 2020, Simulia). Thereby, we offer the flexibility of choosing among many matrix material models and fiber (beam) formulations for the original uncoupled quantities KmandKf. In practice, it is important to note that this problem is inherently non-local. Thus, we enforce the embedding using a “super-element” approach by adding the assembled coupling contributions to the global system. For more detail we encourage the reader to read Steinbrecher et al.’s original work and follow-up work [33, 34, 35]. Once this UEL subroutine is defined, we conduct all network analyses quasi-statically using Abaqus Standard. For convergence, we adopt Abaqus’ default criteria and use automatic stabilization. Here we strictly require the ratio of viscous damping energy (Abaqus output variable ALLSD) to the total strain energy (Abaqus output variable ALLSE) at all steps to be smaller than ≤2%. Similarly, we monitor the coupling constraint violation so that ∥e∥c≤1%at every node, where ∥e∥c=∥gc∥/L∗100% ∈RnfandLis the characteristic dimension of the system. Please see Supplement S2 and Fig. S2 where we conducted a consistency test as a first verification of our implementation. 3 Mechanics of Embedded Fiber Networks Figure 1: Verification Problem. A) An embedded fiber within a solid matrix using both our embedding strategy (fiber-to-matrix coupling) and fully resolved coupling. B) The deformed configuration and displacement field using the two approaches. C) Convergence behavior of the fiber-to-matrix coupling method, where hmatrix is the matrix element size. D) Penalty parameter sensitivity for different fiber-to-matrix element length ratios, where hfiber is the fiber element size. 2.2 Verification Problem To verify our implementation of the mortar-type embedding approach, we consider a single embedded fiber with an end force. We model the matrix as having a length of 5 mwith a square cross-section and an edge length of 1 m. We further assume that the matrix material behaves like a nearly incompressible Neo-Hookean solid with a shear modulus of 50 N/m2. In the center of this matrix, we place a fiber of equal length but with a circular cross-section of radius 0.125 m . For this fiber, we also chose an incompressible Neo-Hookean material model with a shear modulus of 21,730 N /m2 such that the matrix and the fiber have equal bending stiffness. We fix one face of the structure and apply a distributed, shear surface load of TY= 2 N /m2to the free face. We solve this problem using two approaches: a) modeling the fiber as a discrete Timoshenko beam (Abaqus elements B31H and B32H) and embedding it into the matrix (Abaqus element C3D8RH) using the mortar-type embedding approach, and b) modeling both, the fiber and the matrix, using solid elements (C3D8RH), see Fig. 1A. We consider the latter approach as the “ground truth” mechanical behavior. Comparing both approaches, we find that our mortar- type embedding approach achieves a normalized L2-displacement error of ∥e∥<1%, while resulting in significant computational savings (total CPU time 24 seconds versus 2,210 seconds on our workstation with a 36-core CPU at 2.20 GHz). A visualization directly comparing the fiber-to-matrix coupling and fully meshed approach is shown in Fig. 1B. The computational savings stem from the fewer degrees of freedom that are necessary for the embedding approach than for the gold-standard approach. Additionally, we investigate the convergence behavior of the embedding approach with respect to matrix element size hmatrix, see Fig. 1C. Assuming a constant penalty parameter of 100N /m2and a matrix- to-fiber element length ratio of 2.5, we find that the normalized L2-displacement error decreases for decreasing hmatrix until the matrix element size approaches the radius of the fiber element where the error increases again, as previously noted by Steinbrecher et al. [ 33]. Finally, we investigate the sensitivity of the normalized L2-displacement error for various fiber-to-matrix element length ratios and penalty parameters. We find that the fiber-to-matrix element length ratio has little effect on the results of our simulations. Importantly, we do not notice any prominent locking phenomena forhfiber/hmatrix>1.0, even for relatively high penalty parameter values, see Fig. 1D. Please see Supplement S3 and Fig. S3-S5 for additional sensitivity studies of this problem. 4 Mechanics of Embedded Fiber Networks In summary, we affirmed that we implemented Steinbrecher et al.’s approach correctly and found that, consistent with the original publication, it is accurate and highly computationally efficient. 3 A Single Embedded Fiber As a first step toward our goal of understanding the mechanics of embedded fibers and, ultimately, fiber networks, we consider a single embedded helical fiber. For this problem, we chose a solid domain size of 2×1×1 mm and define the helical fiber geometry as x=Rcos (t), y =Rsin (t), z =ct, (9) where R= 0.4 mm andc=l/2πresulting in a single complete loop over the fiber length l= 1.8 mm with t=s/p (R2+c2). Here, sis the arc length. We model the fiber material as linear elastic and describe its mechanical behavior via the Young’s modulus Efand the Poisson’s ratio νf. We chose Ef= 6.5 MPa andνf= 0.495as the default values, representing collagen fibers [ 32]. In contrast to the fiber material, for the matrix material we chose a nearly incompressible neo-Hookean hyperelastic material model as characterized through the matrix shear modulus Gm. We further define the fiber-to-matrix stiffness ratio as α=Gf Gm=Ef 2(1 + νf)Gm, (10) where Gfis the fiber shear modulus. For this single fiber example, we assign a fiber-to-matrix stiffness ratio of α= 5,000. Finally, we subject the embedded fiber to uniaxial extension up to 100% strain. We first examine an uncoupled problem in which we extended the matrix and the fiber in isolation, see Fig. 2A on the left. Next, we model the coupled problem where the helical fiber is embedded in the matrix using the mortar-type approach, see Fig. 2A on the right. When comparing the two cases, we observe that embedding the fiber introduces stress heterogeneity and stress concentrations at the outer boundary of the matrix material, as shown in the stress field. Additionally, embedding fundamentally alters how the fiber is loaded. That is, if we split the fiber’s total strain energy into its bending, stretching, and torsional components, we find that extension of the helical fiber alone initially leads to significant energy being stored in the torsional and bending modes, whereas minimal energy is found in the stretching mode, see Fig. 2B on the left. Only at a strain of 67% does the energy transfer from the torsional and bending modes, into the stretching mode. In contrast, when the fiber is embedded, energy transfers from the torsional and bending modes to the stretching mode even at small strains, see Fig. 2B on the right. At 30% stretching becomes the predominant energy mode. This energy redistribution stems from the matrix constraining the helical fiber’s torsional and bending degrees of freedom, see Fig. 2C for a screenshot of the corresponding interfacial forces between the fiber and matrix. Please see Supplement S4 and Fig. S6 for careful sensitivity studies of this problem and Supplement S5 and Fig. S7 for details on how we compute the strain energies in each loading mode. In summary, embedding even a single helical fiber fundamentally altered its mechanics. That is, embedding alters how the fiber is loaded by redistributing its strain energy from the torsional and bending modes to the stretching mode. Thus, as full fiber networks are assemblies of individual fibers, we anticipate that embedding will have a similar non-trivial effect on the network scale. 4 Embedded Fiber Networks 4.1 Fiber Network Topology & Material Models To generate random fiber networks for our simulations, we use a V oronoi-based approach. To this end, we first generate V oronoi networks from random seeds and then enforce realistic connectivity numbers [ 24]. We do so by removing fiber segments with local connectivity numbers of z≥4until reaching our target values, i.e., ⟨z⟩ ≈3.4for collagen-mimicking networks [14]. Note, the nominal fiber-to-solid volume fraction of our networks follows from ρ=lπR2 f L3, (11) where Rfis the fiber radius, lis the total length of the straight segments, and Lis the cubic domain edge length. Please also note, unless specified otherwise, our networks have edge lengths L= 26 µmand fiber radii Rf= 75 nmrepresenting dimensions that are relevant to collagen fiber networks [ 32,14]. Finally, we include fiber crimp by superimposing sinusoidal undulations along each fiber z- or long-axis, i.e., x=Axsin (ωxz), y =Aysin (ωyz), (12) 5 Mechanics of Embedded Fiber Networks Figure 2: Single Embedded Fiber Example. A) Matrix stress and fiber displacement of the isolated (left) and the embedded fiber cases (right). Here Uis the magnitude of the fiber displacement, Lis the length of solid domain, Gmis the shear modulus of the matrix material, and σxxis the solid stress of the matrix material. B) Fiber strain energy ratios for the isolated (left) and the embedded cases (right). Shear energy is less than 0.4%in both cases and is therefore not included. C) Lagrange multiplier field representing the interface line load between the fiber and the matrix at 30% strain. where z∈[0, li];Ax, Ay=cli;ωx=kxπ/li;ωy=kyπ/li, with kx, ky= [1,2,3, . . .]. Here, the parameter c represents the percentage of the fiber crimp amplitude with respect to the straight segment length li. Our default values arec= 10% ,kx= 2,ky= 1. In those experiments where we study network density, we increase the number of fibers by increasing the number of seeds in our V oronoi-based approach. Here, as in the single helical fiber example, we describe the fiber material as linear elastic via the Young’s modulus Efand the Poisson’s ratio νfand chose Ef= 6.5 MPa andνf= 0.495, again, as the default values. For the matrix, we again chose a nearly incompressible neo-Hookean hyperelastic material model as characterized through the matrix shear modulus Gm. In those experiments where we assume that the matrix material is compressible, we model the matrix via an elastic foam energy function [ 36]. Regardless of our choice of the matrix material, the definition of the 6 Mechanics of Embedded Fiber Networks fiber-to-matrix stiffness ratio remains α=Gf Gm=Ef 2(1 + νf)Gm, (13) where Gfis the fiber shear modulus. Also, where it is convenient to describe the whole network mechanics, we do so through the use of “effective network stiffness” which we define as G=∆σxy ∆γ, Gn=∆σyy ∆γ, E =∆σyy ∆εyy(14) where γandεyyare the shear strain and axial strain under simple shear and uniaxial extension, respectively, and Gnis the so-called normal (Poynting) modulus under simple shear [37]. Figure 3: Embedded Network Mechanics of Varying Densities, Part A. A) Representative networks of rho= 0.1−0.42% density. B) Effect of varying fiber density on average (characteristic) fiber length l/L. C) Shear deformation of a representative embedded fiber network of ρ= 0.42% density at 50% strain. Both the maximum principal matrix stress (ˆσmax) and the fiber axial strain are shown. Here, Gmis the matrix shear modulus. 4.2 Network-Focused Investigation In our first full network example, we model a network of crimped fibers of varying densities and observe their behavior under simple shear. Specifically, we vary the network densities between ρ= 0.0andρ= 0.42%, while setting the fiber-to-matrix stiffness ratio to α= 5,000and the edge length to L= 26 µm, see Fig. 3A,B. Fig. 3C shows an exemplary simulation outcome for a network of ρ= 0.42% density at 50% strain. In our semi- transparent depiction of the simulation domain, it can be seen that fibers align with diagonal direction, i.e., the principal 7 Mechanics of Embedded Fiber Networks direction of strain, but are heterogeneously loaded. The resulting effective stress-strain behavior is shown in both shear and normal direction in Fig. 4A. Toward understanding the mechanics of embedded fiber networks, we first split the total strain energy between the matrix and the fiber network and depict the latter in Fig. 4B. We find that the strain energy contribution of the fiber network increases with network density. We also find that the strain energy contribution of the fiber network increases with strain. When further splitting the network’s strain energy into its bending, stretching, and torsional components we find that the increase of the network contribution with density is driven by increased bending energy, see Fig. 4C. We suspect that this effect is due to reduced characteristic network length that effectively shortens individual fibers and thus increases their resistance to bending. At the same time, we see an increase in energy in the stretching mode with increasing strain, which is likely due to the strain-induced “unfolding” of the crimped fibers as they are strained. Furthermore, Fig. 4D shows the effective shear modulus of the networks as a function of network density and strain. It is clear that the shear modulus increases with strain and that this effect increases with increasing density. In other words, the embedded network induces strain-stiffening behavior in our materials. Similarly, Fig. 4E shows the effective normal modulus of our material. Here, similarly, the normal modulus increases with strain and with network density. Note that by our convention a positive modulus implies that the material contracts in normal direction. Thus, the embedded network amplifies negative Poynting effect in our material. To test whether these behaviors are specific to simple shear, we also repeat the same analysis under uniaxial extension. Interestingly, we find that the behavior is fundamentally the same under uniaxial extension as under simple shear, see Supplement S6 and Fig. S8. Also, please see Supplement S6 and Fig. S9-S11 for sensitivity studies of this problem regarding mesh refinement, penalty parameter, and size effect. In summary, embedding fiber networks lead to strain-stiffening behavior in the combined embedded fiber network and matrix material, which was driven primarily by fiber stretching rather than bending or torsion. Additionally, embedding fiber networks also lead to a more pronounced negative Poynting effect. Finally, we found that this is true under simple shear and uniaxial extension and thus appears deformation mode independent. In the above studies, we have kept fiber crimp and fiber radius constant at ρ= 10 % and R= 75 nm , respectively. Now, we explore their impact on the mechanics of embedded fiber networks with a focus on strain-stiffening and negative Poynting effect, again, under simple shear. To this end, we consider a network with the density ρ= 0.19%, edge length L= 26 µm, and 642total fibers. Next, we vary the fiber crimp ranging from c= 10 % toc= 20 %, see Fig. 5A, and fiber radii between R= 75 nm andR= 100 nm . Fig. 5B illustrates the strain energy phase diagram for these networks under simple shear. Therein, and comparable to our previous observations, we identify two distinct regimes. A bending-dominated regime and a stretching-dominated regime. In the former, the relative bending energy ratio of the fibers is dominant, i.e., greater than all other components of the fiber strain energy. In the latter, the stretching energy component dominates. Thereby, we find that increasing fiber crimp leads to a bending-dominated deformation. Similarly, we find that increasing fiber diameter also leads to bending-dominated deformation. And we again find that increasing strain leads to a transfer of energy from the bending mode to the stretching mode. As a result of this fiber crimp and radius-induced energy transfers, we find that strain-stiffening is larger at smaller crimps and increased fiber radii, see Fig. 5C. We also find that the negative Poynting effect is larger at smaller fiber crimps, but relatively insensitive to fiber radius, see Fig. 5D. In summary, strain stiffening was increased at smaller fiber crimps and with larger fiber diameters. Additionally, negative Poynting effect is also larger at smaller crimps. 4.3 Matrix-Focused Investigations So far, our analyses were focused on the fiber networks. Thus, we next analyze the above simulations with respect to the stress distributions in the solid matrix. Fig. 6A shows an exemplary simulation outcome for a network of ρ= 0.42% density at 50% shear strain. Clearly, the inclusion of a fiber network leads to non-smooth stress fields, i.e., increased heterogeneity. When reducing the volumetric stress fields to probability density functions (PDFs), the distributions flatten with increasing network density, see Fig. 6B. That is, a higher fiber density leads to stresses being more distributed across the problem domain. Finally, we are curious about the impact of fiber networks on matrix compressibility. To this end, we modify the matrix material to that of a compressible neo-Hookean material and vary its Poisson’s ratio between 0.1 to 0.3. Additionally, we vary the network density between 0 and 0.5 %. We then uniaxially extend the networks to 30% uniaxial strain. Fig. 7A shows the average volume change in the matrix material at 30% strain as a function of Poisson’s ratio and network density. We find that all networks allow for significant volume changes as expected from a compressible matrix material. Interestingly, increasing network density leads to smaller volume changes, i.e., the fiber network resists volume changes in the matrix. Illustrating this point, Fig. 7B shows the probability density functions for the volume of individual 8 Mechanics of Embedded Fiber Networks Figure 4: Embedded Network Mechanics of Varying Densities, Part B. A) Shear (left) and normal (right) stress-strain behavior as a function of network density. B) Fiber network strain energy (S.E.) as a fraction of the total strain energy of the embedded network. C) Fiber strain energy components as fractions of the total fiber strain energy. D) The effective shear modulus against shear strain and network density. E) The effective normal modulus against shear strain and network density. Here Gm0andGn m0are the initial (at 0% strain) shear and normal moduli of the matrix, respectively. finite elements for a low-density and a high-density network. We find that increasing network density not only leads to smaller mean volume but also broadens the volume distribution, i.e., again, leads to more heterogeneity. In summary, including fiber networks lead to heterogenization of matrix stresses; an effect that increased with fiber density. Additionally, including fiber networks also resisted volume changes and further lead to heterogeneity in the matrix material. 5 Conclusion In this work, we implemented and used a fiber network embedding scheme to study the mechanics of semi-flexible polymer networks embedded within an amorphous matrix. After verifying our implementation, we explored the impact of embedding on single fibers, fiber networks, and on the embedding matrix. We found that embedding fundamentally alters the mechanics of single fibers and full networks, as well as the embedding matrix. Our investigation focused representatively on two well-known phenomena in the mechanics of semi-flexible polymer networks: strain-stiffening and negative Poynting effect. We found that the fundamental impact of embedding stems from changing the modes in which strain energy is stored in individual fibers and, ergo, fiber networks. That is, while individual, crimped or helical fibers deform largely under bending and torsion, embedding said fibers translates those energies to the stretching mode. The relative contribution depends on both strain as well as geometric and – in the case of full networks – on architectural parameters. That is, with 9 Mechanics of Embedded Fiber Networks Figure 5: Embedded Network Mechanics of Varying Fiber Crimp and Fiber Radius. A) Orthographic projection of the network with 10% fiber crimp (left) and 20% fiber crimp (right). B) Fiber strain energy phase diagram, identifying bending and stretching regimes as functions of fiber radius and fiber crimp. C) The effective shear modulus against fiber crimp, fiber radius, and strain. D) The effective normal modulus against fiber crimp, fiber radius, and strain. Here Gm0andGn m0are the initial (at 0% strain) shear and normal moduli of the matrix, respectively. increasing strain, strain energy is transferred from the bending mode to the stretching mode. In networks, this transfer is delayed with increasing fiber crimp as well as increasing fiber radius that both promote the bending mode. Through careful analysis, we also confirmed that in composites of amorphous matrix material and a semi-flexible polymer network, the latter promotes strain stiffening as well as a pronounced negative Poynting effect. We also explored the effect of embedding fiber networks on the mechanics of the matrix itself. Therein, we specifically found that increasing network density leads to more heterogeneous stress distributions. We additionally found that fiber networks resist volume changes in their host matrix. For soft tissues specifically our findings are highly interesting. Most importantly, they explain the origin of two well-known phenomena, strain-stiffening and negative Poynting effect, which, both, may have important physiological 10 Mechanics of Embedded Fiber Networks Figure 6: Mechanics of the Embedding Matrix Under Varying Network Densities. A) Maximum principal stress ( ˆσmax) distributions in the matrix at ρ= 0% and ρ= 0.42% network densities. B) Probability density function (PDF) of the maximum principal stress for three representative densities. Here Gmis the shear modulus of the matrix material. Figure 7: Mechanics of the Embedding Matrix With Varying Degrees of Compressibility. A) Average volume ratio at 30% uniaxial strain as a function of network density and solid matrix compressibility. B) Probability density function (PDF) of the matrix element volume ratios for a low density (left) and high density (right) embedded fiber network. Here, νmis the matrix Poisson’s ratio. relevance [ 10]. The former is a critical feature that helps prevent soft tissue injury by endowing soft tissues with the ability to progressively resist external forces [ 38,39]. While less clear, the role of latter phenomenon may be related to preventing soft tissues from expanding and enduring compressive stresses within confined space under shear [40,41]. Additionally, our work revealed that semi-flexible polymer networks endow their host matrix with resistance to compression, which may be important in regulating fluid transport in hydrated materials and prevent collapse of arteries that are often contained within soft tissues [ 42]. Finally, our finding that semi-flexible polymer networks lead to heterogeneity in the most material may have important mechanobiological implications [ 43]. That is, fiber network-induced heterogeneity may impact the regulation of cell fate [44]. Through these investigations, we provide fundamental insight into the mechanics of embedded fiber networks and thus advance our understanding of these important materials. However, our work is, of course, not without limitations. Most importantly, we have limited our analysis mostly to a neo-Hookean matrix and linear elastic fiber material. Much remains to be learned about the interactions of other materials and their impact on our findings. Fortunately, our implementation of Steinbrecher et al.’s embedding approach is agnostic to material model choice. That is, any solid material combination available in Abaqus can be used including hyperelastic, viscoelastic, or poroelastic formulations. Thereby, we and others can further our understanding of fiber networks embedded in amorphous matrices for example toward understanding soft tissues and other biological materials. Finally, please note that our work has been focused – in both the motivation and the choices of constitutive, geometric, and architectural parameters – on biological materials. Specifically, we were motivated by collagenous soft tissues 11 Mechanics of Embedded Fiber Networks and have chosen fiber diameter, fiber stiffness, characteristic network size, and matrix stiffness accordingly. However, this framework and our implementation are not limited to collagenous soft tissues. Thus, others who study fibrous composite materials may equally benefit from our work. That is, those interested in other biological materials, or even composites/synthetic materials, may find use in our work. For example, electrospun materials, as used in many industrial, mechanical, chemical, and energy applications, are equally suitable for our type of analysis [ 45,46]. Another potential application exists in the field of tissue engineering. Specifically, the study, modeling, and analysis of local scaffold buckling with respect to various matrix depositions [ 47,48]. Toward this end, we are making our implementation in Abaqus openly available and also describe our approach in much details such that others can adopt our work. To this end, we refer the reader to our GitHub page as listed at the end of this manuscript. Acknowledgments The authors would like to thank Dr. Soham Mane for all the insightful discussions and advice on modeling fiber networks. Funding for this project was provided by the National Science Foundation grants CMMI-1916663, 2046148, and 2127925, and DMR-2105175. Funding was also provided by Office of Naval Research under grant number N00014-23-1-2575 as well as by the National Science Foundation under grant R21HL161832. This support is gratefully acknowledged. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation. Code Availability The interested readers can find the complete pipeline in our GitHub repository: https://github.com/SoftTissueBiomechanicsLab/Embedded_Fibers_UEL