## 2 Materials and Methods
### Theoretical Remarks  
For obtaining the mechanical properties of the samples from the analysis of their nanoindentation data, models arising from the contact mechanics theory were adopted. The Hertz model and the Oliver and Pharr model are the most commonly used models [57]. The latter is suitable for analyzing the responses of hard materials, while the Hertz model can be applied to predict the behaviors of soft materials during the first part of the loading phase, which is related to the elastic region. Considering the spherical tip of the probe  
Figure 3: Nanoindentation curves during an explanatory stress–relaxation test: (**a**) indentation depth vs. time (\(\delta\)-_t_) curve; (**b**) load vs time (\(P\)-_t_) curve; and (**c**) indentation depth vs load (\(\delta\)-\(P\)) curve. and the flat surface of the sample, the Hertz model stated the relation between \(P\) and \(\delta\) as follows:  
\[P=\frac{4}{3}\cdot E_{eff}\cdot R^{\frac{1}{2}}.\delta^{\frac{3}{2}}, \tag{3}\]  
where \(E_{eff}\) is the effective modulus of the contacting bodies (probe and sample), which is defined as follows:  
\[\frac{1}{E_{eff}}=\frac{\left(1-\
u_{sample}^{2}\right)}{E_{sample}}+\frac{ \left(1-\
u_{probe}^{2}\right)}{E_{probe}}, \tag{4}\]  
where \(\
u_{sample}\) and \(\
u_{probe}\) are the Poisson coefficients of the sample and the probe, respectively, and \(E_{sample}\) and \(E_{probe}\) are the Young's moduli of the sample and the probe, respectively. Typically, a probe is made of a much harder material than a sample (\(E_{probe}\gg E_{sample}\)), and thus, Equation (4) can be rewritten as follows:  
\[\frac{1}{E_{eff}}=\frac{\left(1-\
u_{sample}^{2}\right)}{E_{sample}}. \tag{5}\]  
As regards the viscoelastic responses of materials, different constitutive models have been developed to predict them. In this study, the generalized Maxwell model, which is mathematically described by the Prony series, was selected [54; 55; 56]. The analytical formulation describing the load decay during the holding phase as a function of time is reported below:  
\[P(t)=P_{\infty}+\sum_{k=1}^{N}P_{k}\cdot e^{-\frac{t}{\tau_{k}}}\ \ \ \ \ N=2\,\] (6a) where \[P_{(t=\infty)}=P_{\infty}\] is the equilibrium load, which is the value at the end of the relaxation time, and in this case, it was set as \[t=5\ \mathrm{s}\]. The \[P_{k}\] parameters are the material parameters, while the \[\tau_{k}\] parameters are the time parameters, and these represent the 4 parameters of the Prony series to be optimized through the fitting procedure. Alternatively, Equation (6a) can be rewritten as follows:  
\[P(t)=P_{0}\cdot\left(1-\sum_{k=1}^{N}P_{k}\cdot\left(1-e^{-\frac{t}{k}}\right) \right)\ \ \ \ N=2, \tag{6b}\]  
where \(P(t=0)\equiv P_{0}=P_{max}\) is the maximum load, which is the value at the beginning of the holding period (Figure 2b). Using only one Maxwell element (\(N=1\)), the fitting was unsatisfactory (see Appendix A); therefore, a second order Prony series was selected. This behavior was likely related to the principal constituents of the tested biological tissues, i.e., the highly cross-linked collagen and the elastic fibers [58], which were responsible for the different relaxation times. Considering the relationship between the Young's modulus and the shear modulus \(G\) and the conventional assumption of the incompressibility of the tested samples (\(\
u=0.5\)), we determined the following:  
\[P=\frac{4}{3}\cdot E_{eff}\cdot R^{\frac{1}{2}}.\delta^{\frac{3}{2}}, \tag{3}\]0  
Consequently, the link between the effective [PERSON]'s modulus of the sample (\(E_{eff}\)) and the shear modulus of the sample was resolved using \(E_{eff}=4G\).