Multiscale modeling of skeletal muscle to explore its passive mechanical properties and experiments verification

Abbreviations: RMSE: Root mean squared error; RVE: Representative volume element; ECM: Extracellular matrix; AI: Angular integral method

1.   Introduction

Skeletal muscle has complex microscopic structures and multiple compositions. Muscle fibers (also called muscle cells) are the structural unit of skeletal muscle and the basic functional unit of biological movement. The outside of the muscle fiber is wrapped by the endomysium, 5–20 muscle fibers form a fiber bundle, the fiber bundle is wrapped by perimysium ^[1], and a certain number of the muscle bundles are wrapped by the epimysium making up macro skeletal muscle organs. The endomysium, perimysium, and epimysium are connective tissue membranes (containing glial fibers and elastic fibers) ^[2]. Therefore, such hierarchical structure and diverse compositions lead to the simulation of the mechanical behavior of skeletal muscle being a difficult point for scholars.

At the macroscale, scholars have established phenomenological models to simulate the mechanical behavior of skeletal muscle, for example, Pierrat et al. ^[3], Gilchrist et al. ^[4], Böl et al. ^[5], Heidlauf and Röhrle ^[6], Ehret ^[7], Blemker et al. ^[8]. There were also phenomenological models for soft tissues, which could be used to predict the mechanical behavior of skeletal muscles by adjusting model parameters ^[9]. Phenomenological models achieved accurate results. However, the phenomenological models did not consider the influence of skeletal muscle microstructure and compositions on the macro-mechanical properties.

The multiscale method establishes a bridge linking the microscale and macroscale ^[10] and reflects the micro information to the macro mechanical properties. In recent years, the research on multi-scale modeling developed rapidly, both numerical methods and analytical approaches based on transformation field analysis ^[11]. At the microscale, the research was mainly divided into the following three points: research on component characteristics, and research on fiber microstructure, research on fiber distribution. These researches are conducted based on the representative volume element (RVE).

Research on component characteristics. Skeletal muscle is often viewed as a complex of fibers and extracellular matrix (ECM). Bleiler et al. ^[12] used the isotropic Neo-Hookean energy function to characterize the passive mechanical behavior of muscle fibers and assumed that ECM was a transversely isotropic material. Based on the composition of muscle fibers, Böl et al. ^[13] believed that the strain energy equation of skeletal muscle fibers should be characterized by the Holzapfel model ^[14] and isotropic Neo-Hookean energy function is suitable for ECM. The strain energy function of fibers proposed by Kuravi et al. ^[15] was a modified form of the transversely isotropic exponential strain energy function proposed by Holzapfel et al. ^[16]. With the deepening of research, another assumption was proposed, which was that the ECM contained collagen fibers, so that skeletal muscle tissue turned out to be a complex collagen-fiber-reinforcement composite material. It can be concluded that the extracellular matrix should be expressed as an anisotropic material. Spyrou et al. ^[17] split the strain energy of ECM into an isotropic part and a generally anisotropic part to describe the collagen reinforcement. There were some different constitutive models mentioned in the literature ^[2,8].

Research on fiber structure. There were assumptions about the cross-sections of muscle fibers. Through observing skeletal muscle slices stained with hematoxylin and eosin, Sharafi et al. ^[2] proposed that the cross-sections of skeletal muscle fibers and fiber bundles were approximately polygons. Spyrou et al. ^[18] assumed that the cross-sections of skeletal muscle fibers were the regular hexagonals. Subsequently, they chose Voronoi polygons as the fiber cross-sections ^[17]. Kuravi et al. ^[15] used image segmentation and registration methods to build a 3D RVE model of skeletal muscle. Virgilio et al. ^[19] used an agent-based model to randomly generate new muscle fiber geometries. To simplify the geometry structure of RVE, Röhrle et al. ^[20] believed that the polygonal anatomical structure of muscle fibers can be approximated by a circle with a diameter of about 50–80 μm. Bleiler et al. ^[12] pointed out: on average, a circle can best describe the shape of heterogeneity. It should be noted that the volume fraction (${v_f}$) of health muscle fibers has reached 0.95. However, it is difficult to achieve such a high volume fraction using a geometric model with circular cross-sections ^[21]. Fortunately, the regular hexagons, Voronoi polygons, and some irregular polygon shapes created from histological sections can reach 0.95.

Research on the distribution of muscle fibers and collagen fibers. For muscle fibers, the random distribution is realized by randomly generating seeds, and the random sequential addition algorithm was often used ^[17,21,22]. At present, there are two main methods to express the distribution of collagen fibers in constitutive functions, respectively named the generalized structure tensor method and the angular integral method (AI, also referred to as full network model or microsphere model) ^[23]. Based on the AI, Li et al. ^[24,25] introduced a discrete fiber dispersion model based on the triangular discretization of a unit sphere with a finite number of elementary areas. Bleiler et al. ^[12] introduced a concise model of the collagen fibers based on the AI model and applied it to the distribution of collagen fibers in the ECM of skeletal muscle. For the comparison of the generalized structure tensor and AI method, please refer to the literature ^[26]. However, Sharafi et al. ^[2] proved that the distribution of fiber bundles was not completely random and there was a certain hierarchical structure.

After consulting literature, shear deformation plays an important role in the force transmission of skeletal muscle ^[27]. The absence of strain invariant ${I_5}, {I_7}$ in the strain energy functions will lead to unsatisfactory shear simulation results, that is, the shear modulus of the model in different planes are the same. In order to avoid this situation, we build constitutive models containing ${I_5}, {I_7}$ for muscle fibers and ECM respectively ^[28,29,30]. In this manuscript, constitutive models with ${I_5}, {I_7}$ were proposed for capturing the mechanical properties of skeletal muscle under shear deformation. In addition, microscopic observations show that the cross-section of the muscle fiber is oval ^[13]. To establish a model that is closer to the real structure of the muscle fiber, we proposed curved-edge Voronoi polygons as the cross-sections of muscle fibers.

2.   Materials and methods

Skeletal muscle is always assumed to be a transversely isotropic, compressible, and hyperelastic fiber-reinforced composite material ^[6,7]. However, Böl et al. clearly showed that muscle tissue must not be described as a classical fiber-reinforced material ^[31]. Skeletal muscle has a complex collagen-fiber-reinforcement microstructure ^[12]. In other words, in addition to muscle fibers, skeletal muscle also contains collagen fibers which are distributed in ECM. Therefore, we assume that skeletal muscle is composed of three-phase materials and the ECM is laterally anisotropic material ^[32,33,34].

2.1.   Constitutive model

The deformation of materials is represented by the deformation gradient tensor ${\bf{F}}$. The right Cauchy-Green tensor is given by ${\bf{C}} = {{\bf{F}}^T}{\bf{F}}$ ^[35], and the left Cauchy-Green tensor ${\bf{B}} = {\bf{F}}{{\bf{F}}^T}$. We assume that the skeletal muscle fibers are transversely isotropic materials, and regard the ECM containing collagen fibers as an anisotropic material. The strain energy $W$ of skeletal muscle was the superposition of the strain energy ${W_f}$ of muscle fibers and the strain energy ${W_m}$ of ECM:

2.1.1.   Constitutive model of the muscle fiber phase

The muscle fiber is a transversely isotropic material with a preferred direction expressed as a unit vector ${\bf{M}}$ (parallel to Y-axis) in the reference configuration. Smith et al. ^[36] performed a uniaxial tensile test on single muscle fibers, and the Poisson ratio $\mu = {\text{0}}{\text{.4657}} \pm {\text{0}}{\text{.124SD}}$. We assumed that muscle fibers were compressible, transversely isotropic materials.

At the macroscale, scholars have proved that the constitutive model of transversely isotropic materials not only is a function of strain invariants ${I_1}$, ${I_3}$, ${I_4}$ but should include the contribution of strain invariant ${I_5}$ ^[37,38,39]. Based on previous research, we established a constitutive model for muscle fibers,

where ${f_f}(J)$ characterizes volume part, ${W_{f - iso}}$ represents the isotropic part, ${W_{f - ani}}$ represents the anisotropic part, $\kappa _0^f$, $c_1^f$, $c_2^f$, ${k^f}$, $c_3^f$ are the model parameters, and $J$ is the volume ratio representing the volume change. $I_1^*$, ${I_3}$, $I_4^f$, $I_5^f$ are defined as:

The second Piola-Kirchhoff stress tensor can be directly calculated by derivation of the strain energy function ${W_f}$ with the product and chain rules,

where $W_i^f = \frac{{\partial {W_f}}}{{\partial I_i^f}}, (i = 1, 3, 4, 5)$.

The Cauchy stress tensor also called Euler stress tensor is given as:

2.1.2.   Constitutive model of the extracellular matrix phase

The ECM of skeletal muscles plays an important role in lateral force transmission ^[40,41]. The ECM contains a variety of collagen fiber types, such as type I, III, IV, V, VI, VIII, XII, XIII, XIV, XV, XVIII and XIX ^[42]. Adding multiple types of collagen fibers to the skeletal muscle ECM is feasible and not difficult at the calculation level, but it is a difficult task to determine the optimal direction of each type of collagen fibers because the directions of the collagen fibers are not the same in different positions ^[40]. Type I, III collagen form quantitatively the majority of muscle collagens ^[43]. Type I collagen characteristically forms strong parallel fibers and confers tensile strength and rigidity. It exhibits the largest contribution to the stiffness of biological tissues ^[44]. Here, we only consider the influence of type I collagen, so it is assumed that the ECM is a composite material reinforced by type I collagen. The mean directions of collagen fibers are expressed as ${{\bf{M}}_{\text{1}}}, {{\bf{M}}_2}$ distributing in the X–Y plane, and the angle with the X-axis is ${\theta ^m}$ in the reference configuration, as shown in Figure 1. With the deformation along the fiber direction, ${\theta ^m}$ will change accordingly.

The constitutive model of the ECM is expressed as:

where $\kappa _0^m$, $c_1^m$, $c_2^m$, ${k^m}$, $c_3^m$ are the model parameters, and $J$, $I_1^ *$ are calculated as shown in Eq (2.6), the calculation method of $I_4^m$, $I_5^m$, $I_6^m$, $I_7^m$ as follows:

From Eqs (9) and (10), we got:

where $W_i^m = \frac{{\partial {W_m}}}{{\partial I_i^m}}, (i = 1, 3, 4, 5, 6, 7)$.

The constitutive models of the muscle fibers phase and the ECM phase developed in this section have been implemented in the ABAQUS 2016. A particular "user material subroutine" (UMAT) has been developed through the FORTRAN language. Table 1 shows the pseudo-code of calling UMAT and the calculation in UMAT. Please refer to literature ^[45,46] for the calculation of the Jacobian matrix.

2.2.   Definition geometry structure of RVE

The fiber cross-sections designed in RVE are polygons and circles, as introduced in the Introduction. However, due to the geometric conditions of the circular cross-sections, it is difficult to reach a very high value of the volume fraction of fibers. In the literature, the volume fraction with circular cross-sections of fibers is usually set at 0.2–0.5 ^[21] but the volume fraction of skeletal muscle fibers reached 0.95, so that the circular cross-sections do not accurately reflect the volume fraction of skeletal muscle fibers. In a recent study, Böl et al. ^[13] peeled off the muscle fibers of the pig's hind legs and observed the surface of the muscle fibers as a smooth curved surface through microscope observation. In the research of Gillies et al. ^[47], after the NaOH solution dissolving muscle fibers, the connective tissue showed a honeycomb structure in the scanning electron microscope observation image and the cavities enclosed by the connective tissue is not strictly polygons, that is to say, the cavity interfaces are curved surfaces, which further confirms that the muscle fiber interface is not circular.

Therefore, we designed a curved-edge optimization based on the Voronoi polygons, that is, the Voronoi polygons were transformed into irregular ellipses without edges and corners, named as curved-edge Voronoi polygons, as shown in Table 2. The principle was that mark the midpoint of each side of the Voronoi polygons, and use the B-spline curve to connect the middle of the polygons. Finally, the B-spline curves were the boundary lines of fibers.

To study the effect of volume fraction on the macroscopic performance of skeletal muscle, we set the volume fraction to be 0.7, 0.8 and 0.9 for curved-edge Voronoi polygons and Voronoi polygons, and generated the corresponding geometric models, as shown in Table 2.

2.3.   Periodic boundary conditions

Biological tissues are composed of a finite number of RVE by periodic arrangements. Therefore, the application of periodic boundary conditions is necessary to ensure displacement continuity and force continuity. The periodic boundary condition is expressed as ^[48]:

where ${\bf{u}}$ is the microscopic displacement vector, the second-order tensor ${\bf{\bar F}}$ denotes the macroscopic average deformation gradient of RVE, ${\bf{X}}$ denotes the spatial coordinates in the reference configuration, $\delta$ is the second-order identity tensor, ${{\bf{u}}^ * }$ is a periodic displacement field ^[49,50]. In this paper, we use the open-source ABAQUS plugin(easyPBC) to add periodic boundary conditions ^[51].

2.4.   Biomechanical experiment

To obtain more comprehensive biomechanical properties of skeletal muscle, we conducted uniaxial tension and shear experiments on skeletal muscle. The experiments mainly contained the tension along the fiber direction and perpendicular to the fiber direction, longitudinal shear, and in-plane shear tests. The loading directions for experiments are shown in Figure 2. In Figure 2(a)–(c), the red lines represent the fiber whereas and the yellow area represents the extracellular matrix, in Figure 2(d) the red polygons denote the cross-sections of muscle fibers.

2.4.1.   Experiment preparation

The samples were from the right hind limb of a 4 ± 0.5-month-old female pig (n = 1) collected in the slaughterhouse, so it does not involve ethical review and approval. After transporting to the laboratory at 37℃ in a constant temperature box, the samples were cut to the design rectangle and divided into four groups, which were respectively subjected to tension along the fiber direction (G1, n = 5), tension perpendicular to the fiber direction (G2, n = 5), longitudinal shear (G3, n = 5) and in-plane shear (G4, n = 5). The specimen dimensions are shown in Table 3. The data in Table 3 are the average of the measured values of three positions of test areas. Then, these samples were stored in Krebs- The force sensor selected Zwick Roell KAP-TC (500N). To ensure quasi-static stretching, all material tests were performed at a strain rate of $0.5\%/s$ to minimize the viscoelastic effect ^[52,53]. In order to avoid the rigor mortis effect, the whole experiment was done within 7 h post-mortem.

2.4.2.   Uniaxial tension test

In the uniaxial tension tests, the initial distance between the upper and lower clamps is 10 mm as shown in Figure 4, and the beam displacement of the electronic material test device was used as the tension displacement of samples. All uniaxial tension tests were performed up to failure of samples.

The displacement $u$ of the upper clamp was recorded. Nominal strain can be calculated:

where, ${l_0}$ is the initial length, here ${l_0} = 10{\text{mm}}$. The force ${\bf{P}}$ was measured in the tension process and then calculated to obtain the nominal stress ${\bf{S}}$ was calculated as:

where ${A_0}$ is the initial area.

2.4.3.   Shear test

In the skeletal muscle shear experiment, the skeletal muscle samples were installed on a special fixture, as shown in Figure 5(a). The simplified representation of the shear test is shown in Figure 5(b). To exclude the influence of the fixture on the experimental results, we set up a blank control group (n = 5) where all settings were as same as in the shear experiment, but no sample was installed. The measurement results of the blank reference group were averaged, and the measurement results of the G3 and G4 groups were subtracted from the average value as the final shearing experiment result. During the experimental test, if the clamping force is too large, tissue failure will occur at the holding position, and if the clamping force is too small, the sample will slip during deformation. To solve this problem, we laid a layer of 220 grit sandpaper between the fixture and samples.

According to the representation shown in Figure 5(b), the nominal shear strain ${\bf{ \pmb{\mathit{ \pmb{\mathsf{ \gamma}} }} }}$ is calculated:

where $\Delta {\bf{L}}$ is the shear displacement, $n$ is the width of the shear area.

The nominal shear stress can be calculated as:

where, ${\bf{P}}$ is the force applied in the experiment, $L$, $h$ are the length and thickness of the shear deformation area.

2.4.4.   Experimental results

For ease of interpretation, we expressed the tension along the fiber direction as Y tension, and the tension perpendicular to the fiber direction as X tension. Each group (G1–G4) of experiments contained 5 experimental samples, and the measurement data of each experimental sample were smoothed (by Savitzky-Golay). It should be pointed out that our experiment samples were stretched until the sample occurred plastically deformation, and we only researched the elastic deformation range. Refer to the method of Kohn et al. ^[32], the e lastic range is the area where the slope is monotonically increasing. Therefore, in Y tension, we intercept the strain at the range of 0–0.55 as shown in Figure 6. When $\varepsilon > 0.55$, the skeletal muscle tissue was destroyed, plastic deformation occurred and the curve began to fall. Figure 6 shows the average and standard deviation of the experimental measurement data in each group. The nominal stress increases with the increase of the nominal strain under Y tension. When the nominal strain is in the range of 0.3–0.4, the change of the slope of the nominal stress-strain curve increases significantly. When comparing the results of Y tension and X tension, it can be found that the hardness of skeletal muscle along the fiber direction is greater than that of perpendicular to the fiber direction (see Figure 6). This phenomenon also appeared in the experiments of Morrow et al. ^[52]. For shear deformation, the nominal shear stress of the longitudinal shear is greater than the nominal stress of the in-plane share, which is probably caused by the reinforcement of muscle fibers in the longitudinal shear.

3.   Results

3.1.   Definition of the model parameters

The model parameters ($\kappa _0^f$, $c_1^f$, $c_2^f$, ${k^f}$, $c_3^f$, $\kappa _0^m$, $c_1^m$, $c_2^m$, ${k^m}$, $c_3^m$ in Eqs (2)–(5) and (10)–(13) are related to the simulation results of the model and the convergence speed of the model, so it is very important to select a set of reasonable model parameters.

When setting parameter values of muscle fiber in Eqs (2)–(5), we used the strain energy functions of muscle fibers to derive the relationship between nominal stress and stretch ratio. The model parameters were obtained using the tensile experimental data of Böl et al. ^[13] and the volume experimental data of Smith et al. ^[36] fitted with the deduced nominal stress-strain relationship based on the Levenberg‑Marquardt type minimization algorithm ^[54]. For the determination of the model parameters of ECM in Eqs (10)–(13), the same method was used as that of muscle fibers, and the experimental data came from the tensile experimental data of Kohn et al. ^[32] and of Smith et al. ^[36].

When processing the experimental data of Smith et al. ^[36], we used the experimental data fitted with a linear function to obtain the linear relationship of volume ratio-stretch ratio and assumed that the linear relationship passes (0, 1) point to ensure that the volume ratio is 1in the initial configuration. Finally, the functions of ECM muscle fiber were obtained:

where, ${\lambda ^m}$, ${\lambda ^f}$ represent the stretching ratio of ECM and muscle fibers during tensile deformation.

However, the obtained model parameters are referenced values, and the model parameters need to be adjusted according to our experimental data. The final model parameters are listed in Table 4.

Figure 7 shows the comparison between numerical simulation results and experimental data. The RMSE of each deformation has been explained. The RMSE was 0.0035 MPa in the tension along the fiber direction, the RMSE was 0.011254 MPa in the tension perpendicular to the fiber direction, the RMSE was 0.000602 MPa in the in-plane shear, the RMSE was 0.00085 MPa in the longitudinal shear. After analysis, the reason may be that the connective tissue contained in the sample was not completely removed, resulting in high experimental data. The calculation results of the RMSE and Figure 7 prove that our multiscale model can well predict the mechanical response of skeletal muscles during Y tension and shear deformation.

To research the influence of model parameters (sensitivity analysis) on the simulation results, we set $\kappa _0^f$ as 0.0031444, 0.031444, 0.31444, 3.1444, 31.444(MPa), $c_1^f$ as 0.005, 0.006, 0.007, 0.008, 0.009(MPa), $c_2^f$ as 1.4044E-4, 1.6044E-4, 1.8044E-4, 2.0044E-4, 2.2044E-4(MPa), $c_3^f$ as 0.0036, 0.0046, 0.0056, 0.0066, 0.0076, ${k^f}$ as 0.6043, 0.9043, 1.2043, 1.5043, 1.8043. Because the volume ratio of muscle fibers is significantly greater than that of ECM, we only study the influence of model parameters of the fiber phase on the nominal strain.

3.2.   Homogenization response

In In multiscale modeling, periodic boundary conditions satisfy the so-called averaging theorem, which requests that the average of the work decomposes into the product of work-averages ^[55].

where $\sigma$ is the microscopic stress, $\varepsilon$ is the microscopic strain, $V$ is the volume of RVE. To evaluate the effective stress $\bar \sigma$, the stress of the macroscopic level and microscopic level are coupled through the following relation:

where, ${\sigma _i}, \; {\varepsilon _i}, \; {V_i}$ are respectively the stress, strain tensors, and the volume of the ${i^{th}}$ element in RVE, N is the total element number of RVE. Based on the macroscopic Cauchy stress $\bar \sigma$, the macroscopic nominal stress can be obtained $\bar p$ :

Finally, the RVE geometric models were established, the constitutive models were determined, the model parameters were obtained, and the macroscopic stress was calculated. The workflow is shown in Figure 8.

3.3.   Simulation results of uniaxial tension

3.3.1.   Tension along the fiber direction

In this section, we choose the RVE model with ${v_f} = 0.7$ and the Voronoi polygon cross-sections of fibers as an example to study the influence of constitutive model parameters and the simulation results are shown in Figure 9.

When changing $\kappa _0^f$, the nominal stress results are completely consistent except $\kappa _0^f = 0.0031444$ as shown in Figure 9(a). It can be concluded that the influence of $\kappa _0^f$ can be negligible on the simulation results of uniaxial tension along the fiber direction. The influence of $c_1^f$ on nominal stress occurs in the entire strain range. However, when the parameter values differ by 0.001 Mpa, the spacing between the curves is not evident, as shown in Figure 9(b). $c_2^f$ and ${k^f}$ have the same influence area on the nominal stress, the influence area starts from the strain around 0.35, as shown in Figure 9(c), (e). By comparing Figure 9(c), (e), the nominal strain is more sensitive to ${k^f}$. The influence of $c_3^f$ is also in the whole deformation process, and it is most evident in the ankle area, as shown in Figure 9(d). A smaller value of $c_3^f$ will cause a more evident curvature of the ankle area. In summary, it can be seen from Figure 9 that the nominal stress increases with the increase of the parameters $c_1^f$, $c_2^f$, ${k^f}$, $c_3^f$. Therefore, the influence areas of each parameter on the nominal stress-nominal strain curve can be obtained, as shown in Figure 9(f).

In addition to the model parameters, the cross-sections and volume fraction of muscle fibers also affect the stress-strain curve. Subsequently, we select the RVE structure (Table 2) to obtain the nominal stress-nominal strain curves under different fiber cross-sections and volume fractions, as shown in Figure 10. The increase of volume fraction case the increase of the nominal strain, but the effect of volume fraction is not as significant as that of the model parameter. Under the condition of the same volume fraction, the nominal stress of RVE with the curved-edge Voronoi polygons is slightly smaller than that of RVE with the Voronoi polygons. Therefore, it can be considered that the uniaxial tensile test of the fiber direction is sensitive to the fiber cross-sections.

3.3.2.   Tension perpendicular to the fiber direction

In the tension perpendicular to the fiber direction (X-axis), the tensile load is added in the X-axis direction. The collagen fibers in ECM are distributed in the X-Y plane. Collagen fibers resist displacement load during the X-axis tension process and increase the stiffness of the X-axis direction. Simulation results show that only $c_1^f$ has a significant impact on nominal stress, as shown in Figure 11(a), the other parameters have no impact as $c_2^f$, i.e., all the curves are completely coincident, as shown in Figure 11(b).

The influence of $\kappa _0^f$ and $c_1^f$ is mainly reflected in the effect of the isotropic part of the fiber phase. $c_2^f$, $c_3^f$, ${k^f}$ mainly reflect the influence of the anisotropic part of the fiber phase. It can be seen from Figure 11(b) that the anisotropic part of the fiber phase has no contribution to the nominal strain under the tension perpendicular to the fiber direction.

Next, the influence of volume fraction and cross-sections of muscle fibers on the simulation results is researched. It can be seen from Figure 12 that the volume fraction has a negative effect on the nominal strain of the X tension, i.e., when the volume fraction increases, the nominal strain decreases instead. This is because the anisotropic part of the fiber phase has no contribution to the tension perpendicular to the fiber direction, and with the content of ECM gradually reducing, the contribution of the anisotropic part of ECM is reduced responsively. Therefore, the nominal stress decreases as the fiber volume fraction increases. It is worth noting that the nominal stress of RVE with the curved-edge Voronoi polygons is significantly greater than that of RVE with the Voronoi polygons. Comparing the experimental and simulation results in X tension, it is found that the simulation curve is lower than the experimental curve, as shown in Figure 7. Based on the comparison in Figure 12, it can be obtained that the simulation result of RVE with the curved-edge Voronoi polygons is closer to the experimental data than that of the Voronoi polygons. The simulation results show that the shape of the cross-sections of fibers affects the results of the tension perpendicular to the fiber direction.

3.4.   Simulation results of shear deformation

3.4.1.   Longitudinal shear

Under the longitudinal shear, the model parameters $c_1^f$ and $c_3^f$ have evident effects on the nominal shear stress, as shown in Figure 13(a), (c). The change intervals of $c_1^f$ and $c_3^f$ are 0.001 MPa, but the effect of $c_3^f$ is greater than that of $c_1^f$. Figure 13C shows that the change of $c_3^f$ leads to a large change of the nominal stress, which illustrates the effect of strain invariants $I_5^f$ on shear deformation. Simulation results show that $c_2^f$ does not affect the mechanical behavior of the longitudinal shear, as shown in Figure 13(b). The same conclusion was obtained for ${k^f}$.

Under the same strain, the greater the volume fraction, the greater the nominal stress, as shown in Figure 14. However, different fiber cross-sections get almost the same simulation results. When the volume fraction is 0.7, the simulation results of the two kinds of cross-sections are the same. When the volume fraction is 0.9, there is a slight difference. This phenomenon may be caused by the actual volume fraction of RVE with the curved-edge Voronoi polygons being slightly smaller than the calibration value.

3.4.2.   In-plane shear

For the in-plane shear, the direction of the shear force is perpendicular to the fiber direction, at this time, the isotropic part of the fiber and the ECM play a role. It can be seen that the model parameter $c_1^f$ has an impact on the simulation results and this impact becomes more evident when the strain increases in Figure 15(a). From the point of theoretical analysis, the changes of model parameters $c_2^f$, $c_3^f$, ${k^f}$ are meaningless to nominal shear stress. The simulation results proved the above analysis, as shown in Figure 15(b). Here, we only show simulation results with the changing of $c_2^f$.

From Figure 16, we can see that the changes of volume fraction have little effect on the in-plane shear and the influence of fiber cross-section is not evident on the simulation results. However, the influence of the cross-sections is greater than the volume fraction.

4.   Discussion

4.1.   The influence of model parameters and fibers' structure

After the simulation and analysis, we obtained the influence of each parameter on the simulation result, in addition to the influence of volume fraction and cross-sections of fibers on the simulation result. Here we have carried out a comprehensive synthesis as shown in Table 5.

When we determined the model parameters, the maximum value of volume ratio of fiber phase was 1.01518, so that the influence of $\kappa _0^f$ on simulation results was offset by the small volume change. $c_2^f$, ${k^f}$ only affect the mechanical behavior along the fiber direction. The influence of $c_3^f$ on nominal stress is more pronounced in Y tension and longitudinal shear. While the fiber cross-sections affect the mechanical behavior perpendicular to the fiber direction.

The main advantage of the multi-scale framework proposed in this manuscript is predicting the mechanical behavior of skeletal muscle during shear deformation. Here we compare the results of constitutive modes with and without strain invariants ${I_5}, \; {I_7}$, as shown in Figure 17. When $c_3^f = 0$ (without strain invariants ${I_5}, \; {I_7}$), the nominal stress-nominal strain curves of longitudinal shear and in-plane shear completely coincide. The results of our experiments show that there is a difference in the mechanical behavior of skeletal muscle in the longitudinal shear and in-plane shear, as shown in Figure 6. When $c_3^f \ne 0$, the difference is reflected. Therefore, the strain invariants ${I_5}, \; {I_7}$ in the constitutive model have relevant significance.

The volume fraction ${v_f}$ of muscle fibers affects the four types of deformation. In Y tension deformation and longitudinal shear, the skeletal muscle macroscopic stress is positively correlated with volume fraction; in X tension and in-plane shear deformation, the skeletal muscle macroscopic stress is negatively related to volume fraction. Increasing age and skeletal muscle diseases lead to the weakening of skeletal muscle elasticity. One reason for this phenomenon may be the decrease in volume fraction of muscle fibers.

In Figure 12, the results of the RVE with the curved-edge Voronoi polygons are larger than that of the RVE with the Voronoi polygons. In Figure 7, the simulation result of the RVE with the curved-edge Voronoi polygons is smaller than the experimental result in X tension. It can be seen that the result of the RVE with the curved-edge Voronoi polygons is closer to experimental data results and more suitable as the cross-sections of fibers. The curved-edge Voronoi polygons also have certain limitations, i.e., when the volume fraction is above 0.9, there will be a cross between fibers.

4.2.   Comparison of our model with the model of Spyrou et al. ^[17]

In the previous section, we have discussed the influencing factors of the simulation results. Next, we will compare our multiscale model with the model of Spyrou et al. ^[17]. Because Spyrou et al. ^[17] compared their simulation results with the experiments of Morrowa et al. ^[52], we adjusted the model parameters to adapt to the experimental data of Morrowa et al. ^[52]. The model parameters are shown in Table 6. The comparison of the simulation results of the two models were get in Y tension, X tension, and longitudinal shear, as shown in Figure 18.

In Y tension deformation, for $0 < \varepsilon < 0.17$, our simulation results are between experimental results of Morrowa et al. ^[52] and simulation results of Spyrou et al. ^[17]. While our simulation results are smaller than the simulation results of Spyrou et al. ^[17] at $0.17 < \varepsilon < 0.28$. When $0.28 < \varepsilon < 0.3$, the simulation results of our model again coincide with the simulation results of Spyrou et al. ^[17] and experimental results of Morrowa et al. ^[52].

In X tension deformation, our simulation results coincide with the simulation results of Spyrou et al. ^[17] at the range of $0 < \varepsilon < 0.4$; for $0.4 < \varepsilon < 0.78$, our simulation results are between the experimental results of Morrowa et al. ^[52] and simulation results of Spyrou et al. ^[17]. What is certain is that our model is more accurate than the model proposed by Spyrou et al. ^[17] in X tension deformation.

In longitudinal shear deformation, our simulation result is always between the experimental results of Morrowa et al. ^[52] and their simulation results of Spyrou et al. ^[17]. The comparison shows that our model achieves excellent results in shear deformation.

In general, the simulation accuracy of our model in Y tensile deformation is not as good as the model of Spyrou et al. ^[17], but the accuracy of our model has been equal to or higher than the accuracy of their model both in X tension and longitudinal shear deformation.

However, our multiscale model also has a limitation. When generating the RVE model with the curved-edge Voronoi polygons, it is necessary to try repeatedly to obtain the corresponding volume fraction of muscle fibers, which increases the generation period of the RVE geometric model.

5.   Conclusions

In this manuscript, we provided available cross-section (curved-edge Voronoi polygons) of muscle fibers, which was closer to the real geometric structure of muscle fibers; We established constitutive models that can simulate the shear behavior of skeletal muscle from the material property level, instead of the geometric structure of the muscle fibers; We compared the influence of the Voronoi polygons and the curved-edge Voronoi polygons on the simulation results and found that the s cross-sections of muscle fibers affected the behavior of skeletal muscle in X tension and shear deformation, and it had almost no effect on Y tension deformation; We investigated the influence of volume fraction of muscle fibers on the simulation results and gave a reasonable explanation for the elasticity change of skeletal muscle with decrease. In addition, we designed skeletal muscle mechanics experiments to verify the effectiveness of our multiscale model. The comparison results proved that our model can accurately capture the tension and shear behavior of skeletal muscle.

Ethics approval and consent to participate

Ethical review and approval were waived for this study for the following reason. The samples were taken from the slaughterhouse, and the research process did not cause any physical harm to the animals and did not involve cruelty to animals.

Conflict of interest

The authors declare that they have no competing interests.

Acknowledgments

The authors are very grateful for the research supported by the National Natural Science Foundation of China (No.61972117) and the natural science foundation of Heilongjiang Province of China (ZD2019E007)