The models considered are athermal, i.e., thermal fluctuations do not influence fiber mechanics significantly. This is adequate for fibers of such large diameter and persistence length larger than 1 mm, as reported for fibrils [53,54]. Fibers are represented as beams of circular cross section with axial, bending, torsional, and shear stiffness (the Timoshenko model for beams) [55]. Accounting for the bending stiffness of the fibers is important both because collagen fibers have nonzero bending stiffness (directly measured in Ref. [56]) and because the network would be unstable in the initial configuration if fibers would carry only axial forces. The issue of network stability is often disregarded in microscale models of biological tissue. However, according to Maxwell's criterion [57], 3D structures of fibers with only axial stiffness are unstable if the average connectivity number of the network $z\xaf$ is smaller than 6. With $z\xaf$ between 3 and 4 (close to 4), as measured from collagen network images [48,50], and $z\xaf=z=4$ in our models, the structure is subisostatic and has zero stiffness in the unloaded state. The use of subisostatic structures as models for collagen networks is considered acceptable in some works because such structures develop finite stiffness when stretched beyond a critical strain [58,59]. This mimics the experimentally observed toe region of the stress–strain curves, but may lead to errors in kinematics when simulating the heel region. Such errors emerge from the fact that fibers in 3D networks with parameters in the range relevant for collagen deform mainly in the bending mode over the entire physiologically important strain range, as discussed in the literature [26,60].