Bioinspired Fiber Networks With Tunable Mechanical Properties by Additive Manufacturing

Nature has provided perpetual inspiration in generating advanced materials of impressive mechanical functionalities [1]. Among the several fascinating structural materials nature offers, soft fiber networks with random microstructures are abundant in living systems, as they form the primary load bearing framework of the cellular cytoskeleton [2–4], the extracellular matrix [5], and blood clots [6]. These networks are composed of slender fibers that are soft in bending. When subjected to small to intermediate deformations, the fibers primarily bend [7–9] and realign along the direction of maximal principal stress [10–13], causing strain stiffening of the fiber network [14–19]. Interestingly, a residual stress or prestress in the form of uniaxial tension, commonly observed in biopolymer networks [20–23], stiffens them in shear as well [8,24–26], indicating the existence of a complex coupling between tension and shear in these materials. This unique characteristic of strain-induced stiffening makes these networks mechanically adaptive to environmental complexities [27,28], and some attempts were made in mimicking these intricate mechanics with soft synthetic materials [29,30]. The potential use of such bioinspired networks are diverse, for example, in the development of artificial tissue constructs due to the tunability of their mechanical properties to match the properties of organs [29,31] and in the design of flexible electronic and energy devices owing to their adjustable stretchability [32–34].

Light-based additive manufacturing technologies, such as selective laser sintering [35] and digital projection lithography [36], provide ample opportunities to fabricate soft bioinspired networks with highly complex micro-architectures. Although these technologies have been used to fabricate synthetic scaffolds for myocardial tissue [37], cartilage [38], and liver lobule [39], the prior studies were restricted to fabrication of lattice-based models, wherein the structure was spatially periodic. The networks were able to precisely match the number of fibers and connections to the biopolymer networks, but they lacked the intrinsic morphological disorder that biopolymer networks possess. This disorder is useful to achieve desired global mechanics in fabricated networks; for example, tuning the degree of disorder changes the critical strain for strain stiffening [40]. Unlike lattice structures, aperiodic network models cannot be created expediently on the graphical user interface of conventional computer-aided design (CAD) software [41]. Additionally, most of the existing methods to fabricate random fiber networks, such as chemical vapor deposition [42], electrospinning [43], and drop-casting [44], do not replicate a computer model implying limited control over the network architecture. Addressing these challenges is a starting point in truly emulating the intricate geometry and thereby the mechanics of biopolymer networks.

Here, we filled this gap by developing a method to additively manufacture fiber networks from computational models that resembled the structure of type I collagen gels. The models were developed in a scripting environment, meaning that morphological features like the connectivity (i.e., the number of fibers joining at each node), fiber length, and fiber cross section were readily tunable. The numerical model enabled mechanics-based modeling by treating the fibers as beams in finite element software (similarly to Refs. [26,45]), which was used to quantify the small strain shear stiffness under varying levels of uniaxial tensile prestress. We also manufactured a representative fiber network using a commercial three-dimensional (3D) printer, and we built a novel test device to verify that the printed specimen exhibited tension-induced stiffening in shear, similar to both the finite element predictions and networks made of collagen fibers. The experimental results from the tests on the 3D printed specimen were then compared with the results of the finite element simulation.

We developed a workflow to fabricate computer models of fiber networks and developed a test device to demonstrate strain-induced stiffening in a representative 3D printed network. Below, we elaborate on each step.

We created a two-dimensional (2D) model inspired by the structure of type I collagen. Although collagen matrices are 3D, their strain stiffening behavior wherein their tangent modulus depends on applied stress is exactly preserved in 2D network models [7]. We generated Voronoi networks following Refs. [46,47] and introduced further disorder by performing a random deletion of a few fibers, similar to Ref. [26]. In brief, random seeds were scattered within a 2D domain to generate a Voronoi tessellation, with the edges of the resulting polygons being the fibers. The random fiber network thus generated was subjected to further disorder by arbitrarily removing $3%$ of the total number of fibers from its domain (Fig. 1(a)). As expected from the literature [48], the fiber lengths in this Voronoi-based network were Poisson distributed, and they were fully characterized by the average fiber length (L_f), which in the present case was 3.9 mm. In view of the fact that collagen networks exhibit average nodal connectivity well below Maxwell’s isostatic threshold of 2d (with d being the system dimension) [49], the diluted 2D Voronoi networks had an average connectivity of 2.8, which was less than Maxwell’s 2D isostatic threshold of 4.

Network models like this one, having randomness and tunability in fiber length, alignment, and connectivity, are not feasible to develop using commercial CAD software, because such software is restricted to periodic lattice-based models, as in prior studies [37–39]. Here, rather than using CAD software, we created the line diagram of the disordered fiber network (Fig. 1(a)) in a scripting language. Our steps to convert the line diagram into a fabricable model are described later, and our software is freely available at the link given in the Code Availability Statement.

Individual fibers were modeled as three-dimensional beams of rectangular cross section having width w = 0.5 mm in the x–y plane and height h = 1.5 mm in the z direction. The susceptibility to fiber bending is commonly measured by the ratio of the fiber’s bending stiffness to stretching stiffness, which we refer to as the dimensionless bending stiffness, $κ=EI/EALf2$⁠, where E, A, and I are the Young’s modulus, cross-sectional area, and moment of inertia of each fiber [50–53]. Here, for the rectangular cross section, $κ=w2/12Lf2=1.4×10−3$⁠, which is a typical value for fibers in networks of type I collagen and fibrin [7,24,54].

We used Standard Tessellation Language (STL) to create the geometry of each fiber as an isolated three-dimensional triangular mesh meaning that the surface contours of each fiber were defined by a series of connected triangles. The meshes defining each fiber were subsequently placed along the edges of the polygons in the Voronoi network (Figs. 1(b) and 1(c)). Finally, the meshes were combined and exported as a single file that contained one object of disconnected fiber meshes that overlapped one another near the nodes of the network (Fig. 1(c)).

The file from the previous step was imported to the commercial CAD software netfabb (Autodesk) to perform a union of the volumes of the disconnected fiber meshes (Figs. 1(c) and 2(a)) at the regions of their overlap. First, the surfaces of the fiber meshes inside the unified volume of the network were removed. Then the triangles of the disjoint meshes were stitched through re-triangulation to generate a watertight, continuous mesh representing the entire network (Fig. 2(b)). Each edge in this final mesh was shared exactly by two neighboring triangles. The mesh triangles were unique, not self-intersecting, and their normals were always oriented away from the volume enclosed by the mesh. The final mesh, encoded in STL, was ready to be exported to a commercial slicing software specific to the 3D printer’s manufacturer. The necessary scripts and the workflow to generate models of fiber networks and export the STL files for additive manufacturing are available in a public repository at the link given in the Code Availability Statement.

The STL file of the model was processed in a slicing procedure to generate a text file, called G-code, with commands to run a 3D printer for printing the model network. This procedure analyzed the model layer by layer to add to it the necessary mechanical supports during the printing process, which later were manually removed from the finished printed part.

The network was printed by a commercial 3D printing service provider (Midwest Prototyping, Blue Mounds, WI) using an M2 3D printer (Carbon Inc.), which is based on Digital Light Processing (DLP) technology. In brief, the DLP method rapidly fabricates a 3D object layer by layer through spatially controlled solidification of a photo-curable resin [55,56]. Due to high achievable in-plane resolution as low as 75 μm and out-of-plane resolution ≤ 100 μm, this method is suitable for printing complex lattices as demonstrated by prior studies [57,58]. Our goal was to fabricate a soft network (similar to Ref. [59]) with fibers undergoing bending-dominated deformation under loading, and therefore, the resin chosen was a shore A polyurethane elastomer, EPU 40 (Carbon Inc.), with durometer hardness 68 as reported by the manufacturer. The values of Young’s modulus and Poisson’s ratio of EPU 40 have been reported to be 6.81 MPa and 0.48, respectively [60].

The in-plane (x–y) dimensions of the print bed of an M2 3D printer, 189 mm × 118 mm, constrained the design dimensions of the model network to be printed. We chose to print a square model having overall in-plane dimensions 110 mm × 110 mm and out-of-plane thickness 1.5 mm (Fig. 3(a)). The densified regions at the top and bottom of the network (labeled “grip” in Fig. 3(a)) facilitated gripping the printed specimen (Fig. 3(b)) on the test frame during the mechanical testing described in the Experimental Device to Mechanically Test the Printed Specimen section.

We applied simple shear combined with uniaxial prestress on the numerical model of the network by performing 2D finite element simulations using the commercially available finite element software abaqus (Dassault Systèmes). Fibers of rectangular cross section (0.5 mm × 1.5 mm) were modeled as elastic and isotropic beams using two 3-noded quadratic beam elements, following earlier studies [51,61–64]. Some simulations considered simplified linear material response for the fibers (Fig. 4(a)), with the value of Young’s modulus E = 6.81 MPa, as reported previously for EPU 40 [60]. Another set of simulations considered the fibers to possess a nonlinear material stress–strain curve, acquired by digitizing the manufacturer’s published stress–strain curve for EPU 40 (Fig. 4(b)). In all simulations, the Poisson’s ratio (ν) of the fibers was set to 0.48, as reported previously [60]. The connections between fibers were “welded” meaning they transmitted both forces and moments. We evaluated the response of the network using two different solvers as described in the following sections. The global stress and strain measures used in the simulations were the engineering stress and strain.

For instances where the fibers were modeled as linear elastic beams (Fig. 4(a)), the implicit dynamic quasi-static solver with the option of nonlinear geometry was used for the computations, as in our prior work [64–66]. Since the system contained severe geometric nonlinearities, we applied the displacement boundary condition in small load steps (⁠$0.05%$ global strain per step) to ensure that the solver converged to static equilibrium with practically negligible convergence bias. We subjected the network to several levels of uniaxial pretension, $ε=0,10,20,30,40%$ and small simple shear up to $γ=3%$⁠. In terms of the computation time, this solver was ≈5 × faster than the explicit dynamic solver described in the next section.

The tensile stress–strain response of the printing resin (EPU 40) is nonlinear (Fig. 4(b)) and depicts strain softening (Fig. 4(c)). We introduced this material nonlinearity in our model of the network by considering the fibers to possess nonlinear tensile response, exactly matching that of EPU 40 (Fig. 4(b)). The implicit solver described earlier was unable to converge when using the nonlinear material model, so we used the abaqus explicit dynamic solver with the nonlinear geometry option, following our earlier study [13]. Sufficient damping and an optimum time of analysis were considered to attain a quasi-static, steady state at each increment in displacement. Since convergence of the solver can be erroneous if the model is overdampled, we tuned the damping coefficient to ensure that the model was slightly underdamped, meaning that the quasi-static steady state was reached when the model’s kinetic energy (KE) reached a steady state after undergoing a few oscillations. The network was again subjected to several levels of uniaxial pretension, $ε=0,10,20,30,40%$ and small simple shear up to $γ=3%$⁠. At the point of convergence, the model’s steady-state kinetic energy always remained less than 0.001% of its strain energy. See Appendix  A for more details about how convergence of the explicit solver was verified.

The influence of uniaxial pretension on the small strain shear stiffness of fiber networks is well established [8,24,26,64] and is one of the key factors used to assess how well the bioinspired printed network reproduces the desired mechanical response. Here, we developed an in-plane test device to apply uniaxial pretension and simple shear on the printed network (Figs. 5(a) and 5(b)).

The existing in-plane shear testing protocol for a thin specimen, ASTM D 7078 [67,68], cannot induce uniaxial pretension on the specimen necessitating the need for a customized test device. The device, shown conceptually in Fig. 5(a), consisted of two-quarter-inch thick aluminum plates (labeled “plate 1” and “plate 2”) to which the edges of the printed specimen were glued using a cyanoacrylate adhesive (Loctite 401). The dense grip edges of the specimen (Figs. 3(b) and 5(a)) increased the contact area between the glued surfaces and the metal plates and, hence, practically eliminated slip during loading. The use of grips that fix the displacements along with a square specimen can produce a boundary effect, which is common in fiber networks [61,63,69]. We considered reducing the boundary effect by increasing the spacing between plate 1 and plate 2 to use a specimen with an aspect ratio larger than unity. To this end, we simulated specimens having varying aspect ratios, but the results showed that even an aspect ratio of 6 would have a substantial boundary effect, meaning the boundary effect could not be eliminated for any dimensions that would be feasible to print (Appendix  B). For this reason, we chose to use a square network for the remainder of our study.

The y position of plate 2 was adjusted before testing to induce a desired level of uniaxial pretension in the specimen. At a given pretension, the device fully constrained the movement of plate 2 and allowed the plate 1 to slide on a well-lubricated ball bearing plate in the positive x direction. The translation of plate 1 along y was restricted by the tension produced in three long (≈2 m) wire ropes connected to it. All the three wire ropes carried equal tension forces (T₁ = T₂ = T₃) when the specimen was under only pretension. Shear deformation was induced by displacing plate 1 in the positive x direction using a micrometer (Mitutoyo, spindle pitch 0.5 mm) connected to the plate through a load cell (Phidgets Inc., 0–100 g). Given that our objective was to quantify the linear shear stiffness in response to different uniaxial pretensions, the range of shear strain applied was 0 to $≈5%$⁠. If a larger shear strain were desired, the only changes to our test device that would be required would be to use a different micrometer and load cell. Shearing tended to cause rotation of plate 1, but the long wire ropes attached to plate 1 remained nearly parallel to the y axis and supported unequal tension forces, T₁ ≠ T₂ ≠ T₃, to restrict its rotation and maintain static equilibrium. A photograph of the test device is shown in Fig. 5(b).

Calibration of the load cell was necessary before testing. By construction, this load cell was a cantilever beam with four strain gauges attached to it in a Wheatstone bridge arrangement. At a given applied shear displacement on the specimen, the live end of this cantilever beam deflected causing the electrical resistance of the strain gauges to change, unbalancing the Wheatstone bridge. This circuit was energized with an input voltage, and the ratio of output to input voltage was read by a high-resolution analog-to-digital converter (Wheatstone Bridge Phidget, Phidgets Inc.) and recorded by a python script through the use of the Phidget22 software library available on the manufacturer’s website. To calibrate the load cell, we used several combinations of standard calibration weights. The fitted calibration curve relating the force and the voltage ratio was linear up to a force of 2.2 N (Fig. 6), meaning that this load cell performed linearly even at loads 2.24 × the maximum prescribed capacity of 100 g. We used this calibration curve for the tests on the fabricated specimen.

To verify that the fiber network exhibited pretension-induced shear stiffening, we began with the simplified finite element model of the fiber network considering the linear, elastic constitutive model for the fibers (Fig. 4(a)). For this simplified model, we used the implicit dynamic quasi-static solver for its robust convergence and computational efficiency, as discussed in Methods section. We induced uniaxial pretension followed by simple shear on the network. The boundary conditions were applied on the subset of nodes at the top and bottom grips of the network (highlighted blue in Fig. 7(a)). The nodes in the bottom subset were fixed, and the nodes in the top subset were translated in the positive y direction during pretension (⁠$ε$⁠). Subsequently, considering the prestrained configuration as the starting point, simple shear (γ) was applied by translating the nodes in the top subset in the positive x direction (Fig. 7(a)). A representative deformed configuration, corresponding to pretension $ε=30%$ and shear $γ=3%$⁠, is depicted in Fig. 7(b) with the fibers colored by the magnitude of displacement, |U|. We determined the stress–strain (τ–γ) responses in small shear (up to $γ=3%$⁠) at different levels of pretensions $ε=0$⁠, 10, 20, 30 and $40%$ (Fig. 7(c)). These curves were all linear and the slopes, i.e., the values of small strain shear stiffness G₀, increased with increasing pretension. G₀ increased by a factor of 12.1 starting from the value of 82.2 kPa at zero pretension (⁠$ε=0%$⁠) to 998 kPa at the highest pretension (⁠$ε=40%$⁠). These observations were similar to observations of prior studies [8,24,26,64], and they confirmed the occurrence of pretension-induced shear stiffening in the network.

Significant axial tensile strains were present in the fibers at large pretensions of the network. For example, the average axial strains experienced by the fibers were $≈16%$ and $≈23%$ at pretensions (⁠$ε$⁠) of $30%$ and $40%$⁠, respectively (Figs. 7(d) and 7(e)). The resin used to print the specimen (EPU 40) softens at these high levels of strain (Fig. 4(c)), meaning the use of a linear material model for the fibers could cause errors in predictions of the pretension-induced shear stiffening in the specimen. Therefore, we updated the model by matching the nonlinear tensile response of the fibers to that of EPU 40 (Fig. 4(b)) and used the explicit dynamic solver as described in the Methods section. We verified that the explicit solver converged correctly (Appendix  A) and repeated the simulations of Fig. 7(c). We again obtained the stress–strain responses in small shear (up to $γ=3%$⁠) at all pretensions $ε=0$⁠, 10, 20, 30, and $40%$ (Fig. 7(f)). The data were linear, and the slopes increased with increasing pretension, but by a smaller factor of 2.64 starting from the value of 82.2 kPa at zero pretension (⁠$ε=0%$⁠) to 216 kPa at the highest pretension (⁠$ε=40%$⁠). A comparison of the results in Figs. 7(c) and 7(f) indicates that strain softening in the fibers suppressed the level of pretension-induced shear stiffening in the network, but did not eliminate it.

With the simulations predicting a stiffening response in the finite element model, our next step was to ascertain how well the fabricated model inherited this mechanical behavior. We began with the simplest situation, with no applied pretension, by positioning plate 2 in the test device so as to subject the specimen to no pretension. The specimen was gradually loaded in simple shear by displacing plate 1 along the positive x direction. Due to compliance of the test frame and load cell, the actual shear displacement undergone by the specimen was smaller than that applied to plate 1. Given that reliable quantification of the compliance of the test frame and the load cell was challenging, specimen displacements were measured independently. To this end, the test was videographed to obtain sequential digital images of the specimen (e.g., Fig. 8(a)). We tracked the images of the moving grip of the specimen (top red box in Fig. 8(a)), and the displacement of its centroid along x, in pixels, was converted to the physical unit of displacement (mm) to get the applied shear displacement u on the specimen. The shear forces (f) recorded by the load cell were plotted against the corresponding shear displacements (u), and the best-fit response was linear (Fig. 8(b)).

We repeated this simple shear test by subjecting the specimen to a nonzero pretension $ε=10%$ by adjusting the position of plate 2 in the test device before applying shear. Subsequently, this procedure was repeated for additional different levels of pretension $ε=20$⁠, 30, and $40%$⁠. One representative image of the specimen is shown in Fig. 9(a), where the specimen was subjected to a pretension $ε=30%$ and shear $γ=3.6%$⁠. The best-fit force–displacement and the engineering stress–strain responses in small shear were all linear (Fig. 9(b)). The slopes of these curves, that is, the values of small strain shear stiffness (G₀) of the network increased with the increasing pretension. With the range of pretension tested, G₀ increased by 2.65 × starting from the value of 77.0 kPa at zero pretension (⁠$ε=0%$⁠) to 204 kPa at the highest applied pretension (⁠$ε=40%$⁠). These observations confirmed a strain-induced stiffening in the fabricated network.

We compared the values of small strain shear stiffness (G₀) of the fabricated network at different pretensions with the corresponding finite element results (Fig. 10). Firstly, we compared the experimental results with the simulations on the simplified network model where fibers exhibited linear material response (Fig. 10(a)). In the absence of pretension (⁠$ε=0$⁠), the stiffness obtained from the experiment (77.0 kPa) matched closely with the simulation (82.2 kPa), but at finite pretensions (⁠$ε>0$⁠), the fabricated network exhibited far less stiffening than the prediction in the simulations.

Next, we compared the experimental results with the simulations performed on the updated network model wherein the fibers possessed the nonlinear material response of the printing resin (Fig. 10(b)). Interestingly, the stiffness obtained from these simulations matched closely with the experiments, both in the absence and the presence of finite pretensions. For example, the values of shear stiffness predicted by the simulations deviated only by $6.8%$ and $5.9%$ at pretensions of $ε=0$ and $ε=40%$⁠, respectively. Even though the resin itself softened, the stiffening caused by the geometric nonlinearity had a larger effect, thereby producing stiffening of the overall network.

Here, we developed a new workflow to generate STL files of fiber networks for additive manufacturing. This method can generate aperiodic, random fiber networks whose morphological disorder can be tuned to match the micro-architectures of biopolymer networks. To demonstrate the technique, we fabricated one representative network using DLP technology in a commercial 3D printer. Both simulations and experiments demonstrated the presence of tunability of shear stiffness, with shear stiffness increasing with increasing pretension. Strain softening of the printed resin partially, but not fully, reduced the amount of pretension-induced stiffening observed.

The workflow to generate STL files described in this article is general in the sense that the fibers can be arranged to form networks of any architecture, having desired distributions of fiber lengths, connections, and alignment. Moreover, this workflow is also applicable for a network of fibers organized in a 3D domain. Fabrication of a 3D network, however, remains challenging, because 3D disordered networks are not self-supporting, meaning that the out-of-plane fibers do not support the in-plane fibers during fabrication, which requires additional supporting scaffolds. The support materials used in commercial DLP fabrication are made of constituents similar to that of the specimen to be manufactured, and manual removal of these supports remains a challenge, meaning that the use of commercial DLP is restricted to 2D networks. This challenge could likely be addressed by complementing conventional DLP with multi-material projection lithography [70], wherein a second resin is used to make chemically dissolvable supports. Therefore, with systematic integration of recent developments in downstream fabrication, and incorporating our workflow to address the existing challenges in upstream modeling, we anticipate that, in the near future, it will be possible to design bioinspired fiber networks for applications requiring controllable mechanical properties.

The code for making the fiber networks is available². The code for making STL files of fiber networks for additive manufacturing is available³. There exist options to generate STL files from both two- and three-dimensional models of networks.

• The National Science Foundation (Grant No. CMMI-1749400).

The authors have no conflicts of interest to declare.

To verify that the explicit dynamic solver converged correctly, we ensured that the model remained slightly underdamped for all simulations. To this end, we studied the time evolution of the model kinetic energy (Fig. 11(a)). The presence of few initial oscillations in the kinetic energy (for time steps <500, Fig. 11(a)) prior to steady state indicated that the solver was not overdamped. The total number of time steps considered in the analysis was sufficient to ensure that the model kinetic energy fell below $0.001%$ of the model strain energy.

Next, we verified the results of the explicit dynamic solver against the implicit quasi-static solver. To this end, we considered the simplified model of the fiber network where fibers possessed linear stress–strain response. Pretension and shear were applied to the fiber network model, and the model was solved using both solvers. The values of small strain shear stiffness predicted by the explicit dynamic solver at all pretensions closely matched the corresponding predictions of the implicit dynamic quasi-static solver (Fig. 11(b)), indicating correct convergence of the explicit dynamic solver.

Both rigidity and proximity of boundaries can alter the mechanical response of fiber networks [61,63,69]. The experiments used a printed square network that was glued to the plates of the test device before applying tension, meaning that the nodes at the boundaries of the specimen were constrained along both the x and y directions, and the simulations reported in the main text matched these conditions. The choice of using a square specimen combined with constraining displacements in the x direction creates a boundary effect, and, in principle, it could be possible to reduce or eliminate the boundary effect by using a specimen with an aspect ratio >1. To investigate whether a specimen of large aspect ratio would exhibit a reduced boundary effect, we performed additional simulations with and without constraining the lateral (x) component of displacement of the boundary nodes (Figs. 12(a) and 12(b)). The values of small strain shear stiffness of the pretensioned model with and without the lateral constraint, G₁ and G₂, respectively, were compared to quantify the effect of the boundaries. For the square specimen (height-to-width ratio H/W = 1), the ratio G₁/G₂ was 1.65 at the highest pretension of $ε=40%$ (Fig. 12(c)). We next studied models of networks with varying height-to-width ratios H/W (Fig. 12(c)). At $ε=40%$⁠, the ratio G₁/G₂ decreased slowly, starting from a value of 1.65 at H/W = 1 and decreasing to 1.34 at H/W = 6. These results suggested that the effect of boundaries was unavoidable for all practically feasible dimensions of the fiber network.