The uniaxial tensile test is the most widely used method for material characterization (Calvo et al., 2010; Martins et al., 2010; Stecco et al., 2013). It provides stiffness measurements through Young’s modulus, and if the sample undergoes loading and unloading cycles, it also offers insights into viscoelastic properties (Peña et al., 2010). Soft biological tissues such as the arteries, heart, and fascia contain fibers oriented in different directions, forming their internal structure. As a result, their mechanical response varies depending on the loading direction (Guo et al., 2023; Ren et al., 2022; Eng et al., 2014). Biaxial tensile tests are commonly used to evaluate the mechanical anisotropy of these tissues (Takada et al., 2023). However, uniaxial or biaxial tests do not always fully characterize deformation states. In certain cases, tissue behavior cannot be solely described as uniaxial or biaxial, making it necessary to include planar tension tests. For example, Acosta Santamaría et al. (2015) investigated the mechanical behavior of the linea alba in the context of laparotomy closure using planar tension tests. For these reasons, in this work, a multidimensional characterization was conducted using UT, PT, and BxT to replicate the strain states in which the fascia primarily functions.

Soft tissues are usually modeled as composite materials consisting of an isotropic base material reinforced by collagen fibers aligned in two different directions (Peña et al., 2010).

To ensure an accurate reproduction of the fascia’s mechanical response, two material models have been considered (Laita et al., 2024). The first model, based on Holzapfel et al. (2000) and proposed by Pancheri et al. (2014), assumes exponential uncoupled volumetric-deviatoric responses and has been widely used to describe the mechanical behavior of fiber-reinforced soft tissues (Peña et al., 2010; Calvo et al., 2010; Eng et al., 2014). The second model, proposed herein, is a modified exponential invariant-based version of the Costa model (Costa et al., 2001), as introduced by Laita et al. (2024), which considers a coupled response. Within the framework of hyperelasticity, both models assume the tissue is incompressible, undergoes large displacements, and exhibits non-linear anisotropic behavior.

A MATLAB script was developed to analyze the optimal combination of tests and optimize the fitting process. Five types of tests were available for fitting (UT, PT, E1, E2, and E3). The number of tests to combine could be chosen while leaving the rest for prediction, in addition to E4 and E5 biaxial ratios. In this way, combinations of three tests were conducted for both uncoupled and coupled models to study the structural parameters obtained by fitting. The model that provides the best fitting and prediction was chosen to study the combinations with different numbers of tests involved.

Given p , a vector of the q unknown parameters of the SEF, the referred minimization problem can be stated as Equation 12: min p ‖ χ p ‖ 2 2 = min p Σ i = 1 N σ i − σ 1 Ψ 2 , (12) where N is the number of considered points, σ is the stress computed from the experimentally measured force, σ Ψ is the analytical stress, q is the number of parameters of the SEF, and the overlined symbols refer to the mean.

For choosing the proper combination p * , we analyze the R-square error, R 2 , the root mean square error (RMSE), ε , and the relative error e r r * (Destrade et al., 2017) of the fit and predictive processes for all possible combinations, as described in Equation 13: R 2 = 1 − ∑ i = 1 N σ i − σ i Ψ 2 ∑ i = 1 N σ i Ψ − σ Ψ ̄ 2 ε = ∑ i = 1 N σ i − σ i Ψ 2 N − q σ ̄ e r r * = max i σ i − σ Ψ p * σ i . (13)

Following the incompressibility hypothesis ( λ 1 λ 2 λ 3 = 1 ) , the analytical expressions for the non-null Cauchy stress terms obtained from our proposed coupled exponential SEF for the biaxial strain state are described by Equations 14, 15: σ l = − 2   C 0   e C 1 I 1 − 3 + C 2 I 4 − 1 2 + C 3 I 6 − 1 2 C 1 − C 1   λ l 4   λ t 2 + 2   C 2   λ l 4   λ t 2 − 2   C 2   I 4   λ l 4   λ t 2 λ l 2   λ t 2 , (14) σ t = − 2   C 0   e C 1 I 1 − 3 + C 2 I 4 − 1 2 + C 3 I 6 − 1 2 C 1 − C 1   λ l 2   λ t 4 + 2   C 3   λ l 2   λ t 4 − 2   C 3   I 6   λ l 2   λ t 4 λ l 2   λ t 2 . (15)

In the case of a uniaxial strain state, the analytical expressions are given by Equations 16, 17: σ l = − 2   C 0   e C 1 I 1 − 3 + C 2 I 4 − 1 2 + C 3 I 6 − 1 2 C 1 − C 1   λ l 3 + 2   C 2   λ l 3 − 2   C 2   I 4   λ l 3 λ l , (16) σ t = − 2   C 0   e C 1 I 1 − 3 + C 2 I 4 − 1 2 + C 3 I 6 − 1 2 C 1 − C 1   λ t 3 + 2   C 3   λ t 3 − 2   C 3   I 6   λ t 3 λ t . (17)

Finally, for the planar tension strain state, the expressions are given by Equations 18, 19: σ l = − 2   C 0   e C 1 I 1 − 3 + C 2 I 4 − 1 2 + C 3 I 6 − 1 2 C 1 − C 1   λ l 4 + 2   C 2   λ l 4 − 2   C 2   I 4   λ l 4 λ l 2 , (18) σ t = − 2   C 0   e C 1 I 1 − 3 + C 2 I 4 − 1 2 + C 3 I 6 − 1 2 C 1 − C 1   λ t 4 + 2   C 3   λ t 4 − 2   C 3   I 6   λ t 4 λ t 2 . (19)

The longitudinal layer is characterized by a high density of collagen fibers forming fascicles, whereas the transverse layer is thinner, as illustrated in Figure 2A. The results demonstrate that fascia is a highly organized tissue with a clearly defined bilayered structure, as shown in Figure 2B. These layers intersect at an angle of approximately 90°. It can be observed that the transverse layer contains only a single row of collagen fibers, a finding consistent with Pancheri et al. (2014). Figure 2C, stained with Picrosirius Red and observed under polarized light, highlights the nearly 90-degree angle between the layers.

Fascia lata, which surrounds the principal muscles of limbs, works preferentially along one direction, with most of the collagen fibers following this preferred direction, which we denoted as longitudinal; hence, the matrix and fiber transversal direction will play a secondary role in the mechanics and functionality of the fascia. Proof of this is the curves for the uniaxial tests shown in Figure 4. For a stretch of λ = 1.055 , the longitudinal behavior is totally different from transverse, while σ 1 has an average stress value of 3.96 ± 1.15 MPa (mean ± STD), σ 2 only achieves a value of 0.60 ± 0.50 MPa. Following the mechanical behavior that soft tissues usually exhibit, the test begins with an initial zone with no stress increment, and then a strain increment appears (toe region). This is because the unfolding fibers are being stretched; when a value of λ = 1.020 is reached, an exponential increase in stress values is experienced.

The planar tension test uses a large aspect ratio between width and length to measure shear properties. According to Moreira and Nunes (2013), for small deformations, the stress–stretch response for planar tension and simple shear is the same. Nevertheless, a divergence between planar tension and simple shear occurs for stretch values greater than 1.30. As we are far from λ values of 1.30, we consider planar tension valid for measuring simple shear properties. Curves for planar tension shown in Figure 5 describe a mechanical behavior with a longitudinal direction that exhibits greater stiffness in contrast to the transversal direction of fibers, as we observed in the uniaxial test. Longitudinal stress values are 4.77 ± 2.45 MPa (mean ± STD), whereas in the transversal direction, we observed 1.13 ± 0.51 MPa for a λ value of 1.072. A less pronounced non-linear behavior is observed compared to uniaxial curves. Regarding the deviation of the longitudinal curves from the mean, it has been noted that planar tension exhibits greater dispersion.

The mean curves depicted in Figure 6 correspond to the last load cycle at each ratio for biaxial tests. The equibiaxial ratio (1:1) exhibits greater stiffness in both the longitudinal and transverse directions than in uniaxial and planar tension tests. Stretching one fiber family implies an increase in the stiffness of the other. Evidence of this effect is clearly observed by comparing the E1 and E4 ratios: Using the equibiaxial as a reference and considering the described effect of the ratios, a greater stretch in one direction leads to a stiffer curve in the opposite direction than the equivalent curve in the equibiaxial ratio, ensuring the proper performance of the biaxial test. This can be observed in Figure 6F, where the mean curves for each direction and ratio are presented.

Table 1 compiles the mean maximum stress and strain values for each direction and ratio obtained. We include the anisotropy ratio ( η = σ l / σ t ) , defined as the ratio between the longitudinal stress value and the transverse stress value for a specific λ value. In order to compare η across the equibiaxial, uniaxial, and planar strain states, a stretch value of 1.037 has been chosen as a reference.

The η for the uniaxial test exhibits the highest value, of 6.08, followed by the η of the planar tension test, which reaches 3.21 and finally, the equibiaxial, where we found a η value of 2.12. The obtained values for η are reasonable given the characteristics of the different strain states, as the equibiaxial test involves both directions. As observed in Figure 6, increased stretching in one direction results in a stiffer curve in the opposite direction. Evidence of this is that the maximum transversal stress, σ 2 , for λ equal to 1.037 is obtained with the equibiaxial test.

A common point observed in all tests is the significant deviation found in the experiments. Two factors contributing to this could be the extraction area, as regions closer to the tendon or bone may exhibit greater stiffness, and the local mechanical demands the tissue must withstand. If one area supports more stress than another, the fiber density must be higher.

Fitting is used to determine the parameters that define the model. It is based on a minimization problem where successive iterations of the parameters are performed until reaching a minimum in Equation 12. The objective of this step is to compare whether the uncoupled or coupled model is more appropriate based on their fitting and prediction capabilities. Figure 7 represents the average experimental curve for the fifth loading cycle (dashed lines) for each direction and the curves obtained from the fitting (solid lines) through the minimization process. Fitting accounts for the entire range of deformation reached in the different biaxial tests. However, for both the uniaxial and planar tension tests, the maximum values of λ only reach 1.04. Thus, all tests are fitted within the same range of deformation.

We derive the parameters for the fitting process by combining three tests. When the uncoupled model (based on Pancheri et al., 2014) was applied, the optimal combination of test with no constraints in the value of parameters was E1, E3 and PT (Figure 7A) with R fit 2 = 0.964 and R pred 2 = 0.879 . Regarding the structural parameters, the following values were obtained: μ iso = 1470   kPa ,   α = 113.47 ,   μ l = 2442 kPa ,   k l = 167.19 ,   μ t = 0.00   kPa ,   k t = 0.01 . We observe that the parameters associated with the family of transverse fibers are equal to zero, which is not physiologically plausible. If we consider the model as structural, there must be a relationship between the parameter and the tissue’s physiology. On the other hand, for the coupled model proposed in this work, the optimal combination was the ratios E1, E2, and E3 (Figure 7b) with R fit 2 = 0.972 and R pred 2 = 0.878 , the values of the structural parameters were C 0 = 13.88   kPa ,   C 1 = 28.78 ,   C 2 = 124.62 , C 3 = 49.07 . Unlike the uncoupled model, the parameter values in this case align with the expected structural function. The parameter associated with the longitudinal direction is greater than that of the transverse direction, and the latter is greater than that of the matrix. Comparing e r r * (see Equation 13) in both models for the maximum R fit 2 obtained with the uncoupled model (E1, E3, and PT), the coupled model exhibits lower relative errors, especially when the transversal direction is fitted, as shown in Figure 8. The longitudinal direction has a similar relative error in both models along λ , but it is slightly lower in the coupled model. The fitting for the uncoupled model was performed while considering constraints to ensure the physical meaning of the parameters. The relative error indicates that the proposed model achieves better results when different strain states are evaluated for soft-fibered tissues with fiber orientations close to 90°. The graphs show that the model better fits stress values for λ > 1.005 . Note that e r r * = 0 means the model perfectly matches the experimental stress.

Observing the better fitting, improved prediction, and the physiological relevance of the parameters, the coupled model was chosen to study how the combination of tests affects the model’s predictive capability, considering its structural nature.

In this section, the combination of one to five tests is analyzed. It is essential to strike a balance between fitting and prediction. When parameters are obtained based on a single strain state, the fitting error is minimal, but the predictive capability is lost as the parameters become specific to that strain state.

Table 2 summarizes the material parameters and errors for each combination. Structural material parameters exhibit similar values, all within the same order of magnitude, except for the first fitting using only one test. As shown in Table 2, fitting becomes more challenging as the number of tests increases and the strain states become more diverse.

Fitting with two strain states ( R pred 2 = 0.854 ; ε = 0.471 ) implies losing precision when predicting tissue behavior for other strain states compared to fitting with three strain states ( R pred 2 = 0.878 ; ε = 0.371 ). It should be noted that an excessive increase in the number of tests used for fitting does not necessarily result in an improvement in prediction error. While increasing from a single strain state to the combination of two may enhance prediction, fitting with four tests ( R pred 2 = 0.861 ; ε = 0.417 ) does not yield a better prediction than fitting with three. In this sense, fitting with one strain state and with five strain states simultaneously was tested to corroborate the previous idea. Using only a single test, the E1 ratio yielded the best results in terms of the physiological meaning of the parameters and the prediction error R pred 2 that was equal to 0.882 with a ε = 0.406; however, the adjustment error R fit 2 was 0.994 with ε = 0.087. Using five tests, the fitting error R fit 2 is 0.958 with ε = 0.247, and the prediction error worsens with respect to the combination of three tests ( R pred 2 = 0,878 ; ε = 0,371 ) with R pred 2 = 0.862 and ε = 0.406.

Figure 9 illustrates the effect of the number of fitting tests on the errors in fitting and prediction. For each number of tests combined, the optimal prediction has been selected; that is, for the combination of three tests, the R 2 and ε values for the E1, E2, and E3 case are depicted.

Figure 10 depicts the prediction curve for the tests that are not included when fitting with the three strain states (E1, E2, and E3). For low strain values, the prediction curve more accurately follows the real behavior experienced in the test. However, it is also observed that the biaxial ratio E5 proves challenging to predict because it represents a strain state that forces greater stiffness in the softer direction of anisotropy, which contradicts the tests used for fitting.

Throughout this work, a multidimensional characterization has been presented, including three different tests that reproduce a wide range of strain states. The results show that fascia is a highly stiff tissue due to its structure, which consists of layers of collagen fibers spatially oriented in two directions. The highest deformation observed in our tests occurs in the plane tension test, reaching a maximum value of 7.5%. In the other tests, the maximum deformation reached is 5%. These elongation values are consistent with those reported in previous studies (Eng et al., 2014; Pancheri et al., 2014; Ruiz-Alejos et al., 2016). Fascia’s mechanical behavior is characterized by high stiffness, especially when compared to other soft tissues such as the myocardium and arteries. This stiffness allows the fascia to sustain high levels of stress with minimal strain, a characteristic typical of collagenous fibrous tissues like tendons. If tension increases by 8%–10%, it leads to visible tearing of tendon fibers, ultimately resulting in tendon rupture (Wang et al., 2012). The difference in stiffness between the longitudinal and transverse directions is related to the thickness and number of collagen fibers in each direction, which are greater in the longitudinal direction than in the transverse direction, as shown in the histological images in Figure 2A). Similar results were reported by Pancheri et al. (2014).

Eng et al. (2014) obtained 3.5 MPa in biaxial tests for a strain of 4%, while in our study, we measured 3 MPa for the same strain range. Pancheri et al. (2014) reported a maximum strain of 6% in biaxial tests and 8% in uniaxial tests. Regarding maximum stress values in uniaxial tests, they obtained 7 MPa for a strain level of 5.5%, whereas in our study, we reached 4 MPa at the same strain level. Comparing stress in biaxial tests, Pancheri et al. (2014) reported 3 MPa for a 4% strain, which matches our results. Additionally, Ruiz-Alejos et al. (2016) found that deep fascia exhibited a stress of 2.5 MPa at 5.5% elongation in uniaxial tests. Both results are within the same order of magnitude, with the difference accounted for by deviation. As observed by Pancheri et al. (2014), the data illustrate that specimens stretched along the longitudinally oriented fibers exhibit higher stiffness than those stretched in the transverse direction. Despite the different origins of the fascia samples, we observed similar values in sheep fascia lata to those reported by Stecco et al. (2013) for human crural fascia under the same stretch range. As seen in the literature and confirmed by our experimental results across different strain states, fascia exhibits high variability. The stress–strain curves presented in this work show that this deviation is consistent with that reported in other experimental studies.

In this study, we evaluated the accuracy of the model proposed by Pancheri et al. (2014). As they described, a generic angle φ is used despite histological sections showing that collagen fibers form an angle between 80° and 90° (Stecco et al., 2009). We proposed a constitutive model based on a coupled strain energy function, assuming a 90° angle between the anisotropy directions representing the fiber orientations in the tissue. This assumption affects the choice of the constitutive model. Referring back to the formulation in Section 2.2, the unit vectors are defined as M = { 1,0,0 } and N = { 0,1,0 } , which, in turn, impacts the expressions used for stress calculation. In the model by Pancheri et al. (2014), the stress value in one direction does not depend on the other, as seen in the expressions for Ψ 1 , Ψ 4 , and Ψ 6 (Equation 8). Although the model can fit the experimental data (Figure 11A), issues arise with the obtained parameters, as they lack structural meaning. Specifically, the parameters related to the transverse fiber direction are reduced to 0, effectively neglecting one fiber direction. When we impose constraints in the minimization problem to ensure that the transverse parameters remain nonzero and greater than those associated with the matrix, the model is no longer able to fit the experimental data properly, as shown in Figure 11B.

To use an uncoupled constitutive model, it is necessary to not assume that the angle between the anisotropy directions is 90°. Instead, this angle becomes an additional parameter in the problem, defining the unit vectors as M = { c o s ( φ ) , − s i n ( φ ) , 0 } and N = { c o s ( φ ′ ) , s i n ( φ ′ ) , 0 } , where φ represents the fiber angle relative to the 1-axis, and thus φ ′ = 90°- φ . As stated in Pancheri et al. (2014), φ is a phenomenological parameter that they compare to the angle formed by fascia collagen fibers, despite describing a structural strain energy function (SEF). With the unit vectors defined in terms of sine and cosine, the analytical expression for stress calculation in one direction depends on the other. The model we propose in this work effectively fits the experimental data while assuming that the fibers form a 90° angle between them. This is because it incorporates both longitudinal and transverse contributions within the same exponential term, allowing stress in one direction to depend on the other. Even if the angle were treated as a parameter, our model could still accommodate it by incorporating it into the vector definitions, providing flexibility in considering different anisotropy angles.

Considering these aspects, the parameter fitting process was optimized using the coupled model proposed in this work. The main objective is to determine the minimal set of deformation states required for fitting in order to obtain accurate parameters that enable reliable predictions of fascia behavior with the fewest possible experiments.

Regarding the optimal combination of tests among the options studied and listed in Table 2, greater emphasis was placed on minimizing prediction error and reducing the number of test types required, as this directly impacts the number of samples and overall testing effort. As shown in Figure 9, which illustrates the variation of fitting and prediction errors with an increasing number of tests, both R fit 2 and R pred 2 stabilize and remain constant beyond three tests. This indicates that including more than three tests in the fitting process does not enhance prediction accuracy. Additionally, the three necessary tests—biaxial ratios E1, E2, and E3—belong to the same test type, reducing the number of specimens required and the overall testing time by eliminating the need for multiple testing machines.

The aim of a computational model is to enable simulations, making predictability a crucial factor. Our proposed coupled SEF demonstrates excellent predictability with only four parameters, considering that it accounts for three strain states. The material parameters we propose for characterizing fascia and predicting various strain states are listed in Table 3.

Throughout this work, we have emphasized the importance of the obtained parameter values in relation to the structural nature of the model used for fitting. There must be coherence between these values and the structural components they represent. In this regard, it is possible to establish similarities with parameters from other studies. The parameters determined in this study represent a solution to a problem that does not have a unique solution. Therefore, direct comparisons of individual values to establish, for example, a stiffness criterion are not meaningful. Moreover, even if the two models are structural, their defining SEFs may differ. In fact, this work presents an SEF distinct from those proposed by Pancheri et al. (2014) and Ruiz-Alejos et al. (2016). Regardless of the absolute parameter values, a clear pattern emerges: parameters associated with the primary fiber direction are greater than those in the transverse direction. In turn, transverse parameters exceed those related to the isotropic component, which corresponds to the tissue matrix and lacks a mechanical function.

This work has some limitations, one of which is that we tested samples from an animal model rather than human fascia. Although our results are similar to those obtained by Stecco et al. (2013), they cannot be directly extrapolated to the human model. Therefore, the parameters we propose should be used with caution in simulations for human studies.

Regarding the coupled SEF proposed in this study, as discussed by Anssari-Benam et al. (2024), the selection of classical invariants for the isotropic component may be suboptimal if I 2 is excluded, and similarly for the anisotropic component if I 5 and I 7 are not considered. The goal of this study is to develop a model that not only achieves a good fit but also enhances predictive accuracy across different deformation states while maintaining a straightforward formulation. To this end, we have chosen to use models that incorporate a simple exponential function and standard invariants commonly referenced in the literature.

Additionally, our model does not account for viscoelastic properties, which play a significant role in the behavior of soft tissues. The viscoelastic properties of fascia are typically analyzed through stress relaxation and dynamic mechanical analysis (DMA), both of which are widely documented in the literature (Bonifasi-Lista et al., 2005; Prevost et al., 2011; García et al., 2012; Calvo et al., 2014). These properties will be the subject of future studies. The perpendicularity of the fibers is considered; however, soft tissues exhibit fiber dispersion relative to the main direction. The next step to enhance the proposed model would be to incorporate a new parameter for dispersion using techniques such as polarized microscopy (Sáez et al., 2016). The mechanical behavior of soft tissues is governed by their underlying microstructure, particularly the extracellular matrix with embedded collagen fibers. Therefore, studying the micromechanical behavior of individual fibers can provide valuable insights into the macroscopic mechanical response. This approach is commonly used in microstructural models, where the behavior of individual fibers is represented and then homogenized by integrating over the surface of a sphere (Alastrué et al., 2009; Gasser, 2011; Weisbecker et al., 2015; Sáez et al., 2016). This work focuses on the macroscopic response and does not account for the micromechanical behavior of collagen fibers.