## 2 Materials and Methods
### Fibril Hyper-Elastoplastic Model  
In this computational frame, the fibril was modeled via an edited version of generalized Neo-Hookean material, where the general expression of the strain energy function takes the following form:  
\[\psi_{fb}(\overline{I}_{1e},\overline{I}_{4e})=\frac{1}{2}\mu^{fb}(\overline{ I}_{4e})(\overline{I}_{1e}-3) \tag{5}\]  
where the total deformation gradient tensor \(\overline{F}=\overline{F}_{e}\overline{F}_{p}\), where \(e\) and \(p\) stand for elastic and plastic deformation [46; 47], \(\overline{I}_{1e}=tr\) (\(\overline{C}_{e}=\overline{F}_{e}\overline{F}_{e}\overline{F}_{e}\overline{F}_{ e}\)), and \(\overline{I}_{4e}=n_{0}\)\(\overline{C}_{e}\)\(n_{0}\)\({}^{\prime}\), with \(n_{0}\) being the reference direction of the fibril in the initial configuration. The shear modulus \(\mu^{fb}\) is a function of the deformation of the elastic fibril:  
\[\mu^{fb}(\overline{I}_{4e})=\mu_{0}\big{(}tanh\big{(}a_{1}(\overline{I}_{4e}-1 )\big{)}-a_{2}exp\big{(}a_{3}(\overline{I}_{4e}-1_{0})\big{)}\big{)} \tag{6}\]  
where \(\mu_{0}\), \(a_{i}\), and \(I_{0}\) represent the shear modulus, dimensionless parameters, and the secondary stiffening of the fibril, respectively. The hyperbolic form of the strain energy function was beneficial to adapt to the evolution of the fibril stiffness predicted by the  
Figure 1: The collagen fibril coarse-grained model. (**a**) Schematic view of the hierarchical structure of the simulation cell. (**b**) Three-dimensional coarse-grained model of a collagen fibril. (**c**) Divalent and trivalent crosslinks in the coarse-grained model. (**d**) Cross-section of the collagen fibril.  
molecular dynamics model [48; 49; 50; 51; 52]. The strain energy function was then employed to connect the effective stress (\(\Sigma_{eff}\)) of the fibril to the yield condition (\(\Phi\)) as follows.  
\[\Sigma_{eff}=\frac{4}{3}\overline{I}_{4e}\frac{\partial\psi_{fb}}{\partial \overline{I}_{4e}}=\Phi \tag{7}\]  
\(\Phi\) is considered a building block function of crosslink density (\(\beta\)) between the tropocollagen molecules. The plastic strain rate (\(\dot{\xi}\)) in the fibril is driven by the single crystal plasticity model [53; 54; 55]:  
\[\dot{\xi}=\dot{\xi}_{0}\Big{|}\frac{\Sigma_{eff}}{\phi}\Big{|}^{1/k}\mathit{ sig}\Big{(}\Sigma_{eff}\Big{)} \tag{8}\]  
where \(\dot{\xi}_{0}\) is the initial plastic strain rate. The flow resistance of the fibril is given by the following equation [54; 55]:  
\[\dot{\Phi}=\chi\ddot{\xi}\left(1-\frac{\Phi}{\Phi_{s}}\right) \tag{9}\]  
where \(\chi\), \(\Phi\), and \(\Phi_{s}\) are the hardening or softening rate, yield strength, and saturated-flow strength of the fibril, respectively. Following, then, the fibril's plastic velocity gradient takes on a deviatoric shape, as seen below.  
\[\dot{F}_{p}\,\,\bar{F}_{p}^{-1}\,\,=\dot{\xi}\,\,\mathit{dev}(n_{0}\otimes n_ {0}) \tag{10}\]  
Under the constraint of Karush-Kuhn-Tucker loading/unloading conditions, this form was numerically integrated to calculate the uniaxial plastic deformation gradient (\(\overline{F}_{p}\)) and then the uniaxial elastic deformation gradient (\(\overline{F}_{e}\)).