Journal the Mechanical Behavior of Biomedical Materials 3D finite element models from serial section histology of skeletal muscle tissue – The role of micro-architecture on mechanical behaviour

In this contribution we create three-dimensional (3D) finite element models from a series of histological sections of porcine skeletal muscle tissue. Image registration is performed on the stained sections by affinely aligning them using auxiliary markers, followed by image segmentation to determine muscle fibres and the extracellular matrix in each section, with particular regard to the continuity of the fibres through the stack. With this information, 3D virtual tissue samples are reconstructed, discretised, and associated with appropriate non-linear elastic anisotropic material models. While the gross anatomy is directly obtained from the images, the local directions of anisotropy were determined by the use of an analogy with steady state diffusion. The influence of the number of histological sections considered for reconstruction on the numerically simulated mechanical response of the virtual tissue samples is then studied. The results show that muscle tissue is fairly heterogeneous along the fascicles, and that transverse isotropy is inadequate in describing their material symmetry at the typical length scale of a fascicle. Numerical simulations of different load cases suggest that ignoring the undulations of fibres and their non-uniform cross-sections only moderately affects the passive response of the tissue in tensile and compressive modes, but can become crucial when predicting the response to generic loads and activation.


Introduction
Skeletal muscles are responsible for generation of force, motion, and gait stability (Barrett et al., 2016;Lieber, 2002). Compared to other soft biological tissues, they can be activated under voluntary control. Structurally, they represent a highly hierarchical composite material system consisting of a multitude of long, multi-nucleated muscle fibres packed into fascicles which, in turn are packed to form muscles (Barrett et al., 2016;Lieber, 2002). Muscle fibres, fascicles, and the whole muscle, respectively, are embedded into a continuous network of collagenous sheaths arranged into endomysium, perimysium, and epimysium that collectively, together with other constituents, form the muscle's extracellular matrix (ECM) (Lieber, 2002). The presence of such intricate micro-structure naturally suggests that it could play a pivotal role in the material symmetry, mechanical response, and physiological load transfer mechanisms in skeletal muscle tissue. However, the precise pathways and the role of the individual tissue components still remain elusive in many regards. For instance, experimental evidence points to a pronounced structural influence on the mechanical response of R. Kuravi et al.  the tissue constituents have dominated over several decades (Lemos et al., 2004;Johansson et al., 2000;Röhrle and Pullan, 2007;Sharifimajd and Stålhand, 2013;Böl et al., 2011;Calvo et al., 2010;Latorre et al., 2018;Seydewitz et al., 2019). More recently, however, approaches which incorporate the particular tissue micro-structure in varied level of detail, have gained increased interest, including both refined continuum models Bleiler et al., 2019) and finite element (FE) based microstructural models that account for several length scales (Spyrou et al., 2017(Spyrou et al., , 2019Marcucci et al., 2019;Virgilio et al., 2015;Blemker, 2011, 2010). A simplistic representation of the muscle tissue, motivated from physiology, suggests a uni-directionally reinforced material consisting of ECM and muscle fibres, which clearly suggests a transversely isotropic material system. This, in fact, is the prevalent material symmetry appearing in numerous continuum mechanical models of skeletal muscle tissue Röhrle and Pullan, 2007;Ehret et al., 2011;Böl et al., 2011;Bleiler et al., 2019;Seydewitz et al., 2019). On the other hand, in micro-structural FE models the non-circular shape and non-uniform distribution of muscle fibre cross-sections can, in general, generate more complex material symmetry (Spyrou et al., 2017(Spyrou et al., , 2019Marcucci et al., 2019;Virgilio et al., 2015;Sharafi and Blemker, 2011) and, in fact, deviations from transverse isotropy have been discussed (Sharafi and Blemker, 2010). A closer inspection of the muscle tissue reveals that (i) muscle fibre cross-sections are typically neither circular nor regular (Briguet et al., 2004;Sertel et al., 2011), and (ii) along the longitudinal direction, fibres have neither uniform shape nor do they follow a straight line (Farrell and Fedde, 1969;Johnson and Beattie, 1973). Moreover, at the length scale of about 100 μm, perimysium winds around the fibres in layers, groups them into fascicles Oshima et al., 2007), and adds even more complexity. The proposed micro-structural FE models of muscle tissue that were used so far to study the effects of the structural and material properties of the tissue components, were either based on artificial cross-sections or single histological sections (Spyrou et al., 2019(Spyrou et al., , 2017Marcucci et al., 2017Marcucci et al., , 2019Blemker, 2011, 2010). Clearly, this disregards any change of micro-structure along the fibre direction. Therefore, the present study is dedicated to investigating the role and the influence of resolving the 3D micro-structure on the mechanical response of skeletal muscle tissue through detailed FE models. To this end, a procedure is proposed to generate fully 3D FE models from a stack of histological cross-section images of muscle tissue. Furthermore, a numerical method is presented to estimate local directions of anisotropy in the collagenous ECM and muscle fibres. These FE models are then aptly parametrised and numerical simulations are performed to study the response of the virtual tissue models to different mechanical loads. The influence resolving the 3D micro-structure is assessed by comparing the response of models generated from 1, 2, and 5 images.

Histology
In line with previous studies of our groups (Böl et al., 2014(Böl et al., , 2016Garcés-Schröder et al., 2018;Böl et al., 2019) histological investigations were performed on porcine biceps femoris muscle. A hind leg ( = 1) of a female domestic pig (Sus scrofa domestica) was obtained from a slaughterhouse immediately after animal sacrifice, and sample preparation started approx. 1 h later, directly after transport of the leg to the lab. Within the next hour, the muscle was excised and cubic samples ( = 5) with a characteristic edge length of 10 mm were cut before being stored for a short time (1-5 min) in physiological saline-soaked cloths. The samples were subsequently stored for 30 min in a climatic chamber (Memmert AtmoControl, Germany) at 4 • C, positioned with an embedding medium on a cork plate and quickly frozen for 20 s in isopentane (-150 • C) cooled with nitrogen (-196 • C). Before further processing, the samples still positioned on the plate were stored in an ultra-low temperature freezer (New Brunswick ™ , Innova ® U101, Eppendorf, Germany) at −80 • C. Before cutting thin sections, three pin holes outside the region of interest (see Figs. 1a,b), were punched into the frozen sample (see Section 2.2.1). Thereafter, thin sections of thickness 14 μm were sequentially cut from the sample using cryo-microtome technology (Leica CM1100, Leica, Germany), placed on microscope slides and stored at room temperature for 24 h. The sections were then stained with Picrosirius red (Appendix A) and digitised using a digital microscope (ZEISS Smartzoom 5) at an image resolution of 5970 × 6632 pixels.

Identification of muscle tissue constituents
The histological information was used to reconstruct 3D geometrical models of the muscle tissue based on the following two main simplifying assumptions: (i) Permanent tissue deformations introduced during the cutting procedure lead to 'artificial' differences between tissue sections. These deformations were approximated and corrected by a global (i.e. homogeneous) affine transformation between the first (reference) and consecutive images, which can be determined from the position of the three pin holes present on all sections (see Section 2.2.1). All differences between the sections that remain after this transformation are R. Kuravi et al. due to the inherent changes of muscle micro-structure along the muscle. (ii) Based on the staining, only two sub-components were identified in each section, i.e., muscle fibres and the matrix surrounding them. Accordingly, the matrix was considered as ECM, and no distinction was made between endomysium and perimysium during segmentation. Likewise, inclusions of the fatty tissue between the muscle fibres were not separately accounted for. To avoid artefacts from an automatised tissue reconstruction, 3D geometries were manually reconstructed from segmented (two-dimensional) histological slices. To create virtual volumes with varying level of detail, subsets of histological sections with = 1, 2, 5 slices were considered.

Image registration using an affine transformation
The image stack of histological sections  = { 1 , 2 , … , } was converted from RGB to binary format using standard plug-ins in Fiji (Schindelin et al., 2012) to accentuate auxiliary markers (pin holes). All three pin holes were then isolated and their respective centroid locations were determined in each image by utilising standard filtering and centroid detection features in MATLAB ® (MATLAB, 2017). These centroids are denoted by O, A, and B in counter-clockwise fashion (Fig. 1a). Vectors and , respectively, pointing from O to A , and from O to B , are then defined on each section = 1, 2, … , of the stack (Fig. 1e, f). Using the first image ( = 1) as a reference, all other images ( = 2, 3, … , ) were shifted by a constant displacement such that the centres of all pin holes O coincided. The four components of an in-plane affine deformation tensor that maps and ( = 2, 3..., ) to 1 and 1 were determined by solving the linear system of equations for each section of the stack as illustrated in Fig. 1. The respective affine mappings af f ∶ → + , where ∈ , defined through and , were then applied to the original RGB images such that to obtain a set of registered images

Identification of muscle fibres and the ECM
Numerous distinct image analysis methods, both semi-automatic and automatic, have been proposed in the literature for muscle fibre segmentation (Wang, 2016;Smith and Barton, 2014), including active contour (Klemenčič et al., 1998;Kim et al., 2007) and ridge-detection based methods (Mula et al., 2012;Sertel et al., 2011), that partly pose requirements on staining and microscopy protocols. In this work, the Xlib plug-in library (Münch, 2019) embedded in Fiji (Schindelin et al., 2012) was used to perform image analysis on the sections stained with Picrosirius red. Muscle fibres and the corresponding ECM portions in each image of the image stack  ′ were gleaned out by subjecting them through a sequence of image processing operations: clustering, thresholding, and disconnect particles in Xlib. A smaller sub-image stack (1500 μm × 1500 μm) selected from the same location of each image in  ′ was cropped and considered for further processing (Fig. 2a). These sub-images, denoted by = {̂1,̂2, … .,̂}, were subjected to clustering based on a K-means parallel algorithm (Kanungo et al., 2002). Clustering effectively segregates muscle fibres and collagen from each other and from other intramuscular components including fat (Fig. 2b). The outcome was then converted to 8-bit format and subjected to a thresholding operation which segregates muscle fibres and all other components into respective white and black areas of the image (Fig. 2c). Noise manifesting as tiny particles in the thresholded images was filtered out utilising the remove outliers plug-in of Fiji. Since the outcome suffers from certain muscle fibres appearing as conjoint, the disconnect particle algorithm (Münch et al., 2006) was applied, that identifies these conjoint muscle fibres and disconnects them at their bottle-neck locations. Finally, each muscle fibre was accorded with a unique colour identifier that distinguishes them from each other (Fig. 2d). The resulting stack is designated as . The parameters chosen for the above operations are summarised in Table 1, whereas the parameters for thresholding operation were chosen according to the image histogram specific to each image.
Although all the above operations could be performed on all the images of the stack , in what follows, five imageŝ ( = 1, 4, 7, 10, 13), equally spaced from each other were chosen for further analyses. They shall be referred to as the image stack and their segmented outcome ( = 1, 4, 7, 10, 13) forms the stack. Table 3 lists the number of R. Kuravi et al.  fibres in each image of the stack and the fraction of fibres that are identified in each image of with 1 chosen as reference image. Although this approach adequately captures all muscle fibres in each image, their cumulative area fraction in each image (0.63 to 0.76, Table 2) is markedly lower than the values > 0.9 reported in literature (Kjaer, 2004). For instance, Bendall (1967) reports collagen content < 10% in various bovine muscles and similar values are reported by Listrat et al. (1999). Lawson and Purslow (2001) report collagen content to be < 1% in the muscles of adult chicken. Notably, the determined muscle fibre area fraction is sensitive to the parameters used for image analysis, in particular for clustering and thresholding. To this end, a further analysis step is proposed, based on a combination of Voronoi tessellation of and low level clustering of.

Adjustment of muscle fibre area fraction
The resemblance between muscle fibres and Voronoi cells was first studied by Honda (1978Honda ( , 1983, and it was used in conjunction with active contour models towards automated segmentation of muscle fibres from histological cross-sections (Klemenčič et al., 1998). However, the initial seeds for Voronoi tessellation had to be provided by an operator by identifying approximate centroids of all fibres manually. Drawing parallels between Voronoi tessellation and muscle fibre disposition in histological cross-sections, Spyrou et al. (2019) numerically generated a statistical representation of the muscle tissue cross-section. The present approach differs from others, e.g. Spyrou et al. (2019) and Klemenčič et al. (1998), in that the centroids of all muscle fibres are readily  Fig. 3c, and Table 2 summarises the corresponding cumulative muscle fibre area fraction, taking values over 90%. The area of the thin lines separating individual Voronoi cells represents the minimum achievable collagen area measured in terms of the number of enclosed pixels, and is hence influenced by the corresponding image resolution. Therefore, at a given image resolution, Voronoi tessellation provides a lower bound for the collagen area fraction in a given image. However, it is observed that thicker endomysium and perimysium bands are concomitantly trimmed by Voronoi tessellation (Fig. 3c). To address this, thick bands of collagen were separately isolated by clustering the image stack with a fewer number of clusters (Table 1), followed by thresholding (Fig. 3d). Finally, a binary AND operation was performed on the outcomes of Voronoi tessellation and low level clustering to adjust the muscle fibre area fraction while maintaining the fibre count (Figs. 2e, 3e and Table 2).

Matching muscle fibre cross-sections
The particle disconnection operation when performed on the outcome of thresholding (Fig. 2c) identifies all muscle fibres with a unique R. Kuravi et al. Fig. 4. Generation of virtual 3D muscle tissue from the boundary coordinates obtained from the segmented stack. (a) All muscle fibre boundaries obtained from images 1 , 7 and 13 . (b, c) Boundary wire splines and corresponding virtual volume of a single muscle fibre. (d-i) All virtual muscle fibres obtained by extrusion, i.e., from 1 (d), from 2 images, i.e., 1 and 13 (e), 3 images, i.e., 1 , 7 and 13 (f), and all 5 images (g). Arrows point at intra-fascicularly terminating fibres. (h) Cubic muscle fibre instance sectioned from (g). (i) The ECM generated from (g). (j) Virtual muscle sample obtained from (h) and (i). colour identifier (Fig. 2d). For a given fibre this colour identifier varies across the stack. Therefore, to track a muscle fibre across the stack, all its cross-sections in were manually mapped before assigning the same colour identifier. To this end, at first, fascicles (i.e., large characteristic features) were mapped across segmented images, and then individual fibre to fibre matches within each fascicle were identified across all images. In difficult cases, auxiliary use of additional images of the stack was made to track the course of single fibres. The fraction of those muscle fibres identified in every image of, termed the fibre continuity fraction, is listed in Table 3. The undulation of a single fibre with respect to the stacking axis of is evaluated by connecting the centroids of its cross-sections across the stack. Column 4 of Table 3 lists the average inclination of all the muscle fibres extending between 1 and ( = 4, 7, 10, 13) with respect to the stacking axis of.

Tissue reconstruction from variable number of sections
By tracking the unique colour identifier of each muscle fibre, the boundary coordinates of its cross-section in each image of were extracted by utilising custom scripts in MATLAB ® (The MathWorks Inc., Natick, MA, USA). The coordinates of 25 equidistant points on each such boundary were introduced as smooth closed wire spline features into Abaqus/CAE 6.14-1 (Abaqus/Standard 6.14-1, 2014a). 3D virtual muscle fibres were then generated from these wire splines by using the create solid loft feature. This process was automatised by custom built scripts implemented in Python 2.7.16. This process is pictorially illustrated in Figs. 4a-c. The start and end tangencies were set to neardefault values of Abaqus so that the loft features approach linearly between first and second and from next-to-last to last section, and hence do not alter the shape or the direction of the loft. 3D virtual muscle fibre assemblies generated from 1, 2, 3 and 5 images of are respectively depicted in images Figs. 4d-g.
Cubic muscle fibre samples were sectioned from the resultant muscle fibre assemblies. Fig. 4h depicts the sample obtained from 5 images. Cubic ECM samples were obtained by subtracting the muscle fibre sample from a cube of identical edge length ( Fig. 4i), and cubic muscle tissue samples were then generated by combining the ECM and the muscle fibre samples to finally obtain the domain of interest (DOI) (Fig. 4j).

Diffusion-based estimation of local anisotropy directions
The aforementioned tissue reconstruction yields a cubic domain encompassing the ECM and muscle fibres, which are anisotropic materials in themselves. In the following, the structural considerations made towards modelling of these materials and the steps involved in deducing their (local) preferential directions are described.

Muscle fibres
Muscle fibres are multi-nucleated cells and exhibit transverse isotropy due to uniaxial alignment of myofibrils. Hence, the axis of transverse isotropy is given by the local muscle fibre direction Blemker, 2010, 2011;Ehret et al., 2011), depicted in Fig. 5a.

Collagen architecture
In skeletal muscle tissues, collagen is the main constituent of the ECM which is embedded in a bath of proteoglycans (Purslow, 2008). The ECM takes the shape of a network-like scaffolding that surrounds muscle fibres and fascicles in a honeycomb-like structure (Trotter and Purslow, 1992;Oshima et al., 2007). It differs both in composition and in alignment between endomysium, perimysium, and  epimysium (Gillies and Lieber, 2011;Purslow, 2008). Perimysium arranges itself around muscle fascicles in multiple layers (Kjaer, 2004). Purslow (1989Purslow ( , 2008 studied the variation of collagen fibril orientation and its waviness in perimysium in bovine sternomandibularis muscle when stretched, and observed that the orientation has two preferred directions at about +55 • with respect to the muscle fibre direction in the rest state. The structure of endomysial connective tissue is less clear. For instance, in cat biceps femoris it was observed to have a relatively homogeneous distribution with no preferential directions which is maintained upon loading (Trotter and Purslow, 1992;Purslow, 2008). However, in mouse extensor digitorum longus muscles, the endomysial tissue was reported to exhibit a longitudinal alignment to muscle fibres at their interface (Gillies and Lieber, 2011). Epimysial connective tissue is layered as well (Gao et al., 2008) where proximity to the muscle tissue dictates the arrangement of collagen network from being randomly distributed to being highly aligned.
In the current work, all collagenous connective tissues were modelled as one type of ECM with properties inspired by the structure of perimysium for the following reasons: (i) The length scale of muscle tissue samples considered is of a few hundred microns which does not encompass epimysium. (ii) The Picrosirius red staining for collagen in histological sections does not allow adequate differentiation between endo-and perimysium connective tissue during image segmentation. (iii) Light et al. (1985) observed that in six different bovine skeletal muscles, the perimysium connective tissue is far more in abundance than the endomysial connective tissue. Therefore, following (Purslow, 1989(Purslow, , 2008, the collagenous matrix is modelled as an orthotropic membrane described by two local in-plane preferred fibre directions 1 and 2 , oriented at an angle + with respect to local muscle fibre direction ′ in a membrane with normal , as shown in Fig. 5b. These vectors in turn can be represented in an orthonormal basis { , ′ , }, where = × ′ . To this end, for a given , it is sufficient to determine the vectors and ′ to completely capture the assumed local orthotropy of the ECM (see also Section 4.3.1). It is noted that the vector ′ is defined locally in the ECM, while is defined locally in muscle fibres (Fig. 5a). They are identical only at the muscle fibre-ECM interface in general, as will be discussed in detail in Section 3.2.2.

Steady-state diffusion
The solution to the Laplacian equation in a region of interest with suitable boundary conditions has been recently utilised to extract fascicle trajectories in skeletal muscles (Choi and Blemker, 2013;Handsfield et al., 2017), and was described as a 3D computational alternative to diffusion tensor magnetic resonance imaging for identifying tissue scale architecture (Choi and Blemker, 2013). A similar approach was implemented in 2D to estimate local collagen fibre orientation in a human atherosclerotic artery (Raina and Miehe, 2016). Appealing to these ideas, steady state diffusion is used in the present work to deduce the local directions of anisotropy in the DOI. To this end, let denote a conservative vector field defined in a domain , and the scalar field variable associated with it related through the diffusion-like constitutive equation where represents the diffusivity of the material. Assuming to be constant, div = 0 yields the Laplace equation with boundary conditions where and , respectively, denote the boundaries associated with Neumann and Dirichlet boundary conditions and represents the normal to the boundary. In order to employ the steady-state diffusion problem (Eq. (4)) to estimate the directions of local anisotropy in muscle fibres (MF) and the ECM, the two tissue components are furnished with two auxiliary (and physically meaningless) diffusion coefficients MF and ECM . For later use, we separate the boundary  of region  occupied by the cubic muscle sample into top (t), bottom (b) and lateral (l) faces. In addition, muscle fibre and ECM portions of the boundary and of each of these faces are also identified resulting in Moreover, all the internal interfaces between muscle fibres and ECM, i.e., all lateral surfaces of the muscle fibres are denoted by  IF , and their corresponding local outward normal vector is denoted by IF . These surfaces are illustrated in Fig. 6. Without loss of generality, steady state heat conduction was chosen for implementation of the above in a FE analysis framework as will be summarised in Section 5.

Muscle fibres
To determine the local muscle fibre directions in the DOI, we consider the problem of steady state diffusion between the top and bottom surfaces of the fibres (Fig. 6b). The boundary conditions are defined as = 1 and = 2 < 1 , with insulated lateral surfaces represented by ⋅ IF = 0 on  IF . The latter condition is approached by a choice of numerical parameters such that MF ∕ ECM ≫ 1, which effectively enforces minimal flux across the fibre boundaries due to relatively low diffusivity of the ECM. By this means, boundary conditions need to be defined only on the outer boundaries of the DOI as where  denotes the outward normal to respective surfaces. The corresponding solution of the boundary value problem (BVP) yields the flux vector , from which the local directions of the muscle fibres are computed.

ECM layers
To estimate the disposition of the ECM layers, we explore its analogy with the concentration field that would result if a constant flux was applied at the fibre-ECM interfaces (Fig. 6e) while maintaining the ECM portions on the lateral surfaces at a constant = 0 . It is assumed that the ECM layers are reasonably well-represented with resulting iso-surfaces of the concentration field , whose local surface normals are given by the flux vector . Hence, a source term IF is introduced at the fibre-ECM interfaces  IF , and by choosing the ratio MF ∕ ECM ≫ 1 it is ensured that each fibre is practically at constant temperature in equilibrium state. The problem is completed by the boundary conditions From the resulting solution of the BVP, local unit normal vectors to the iso-surfaces (ECM sheets) are computed from the flux vector .
In the reference configuration, let denote the local muscle fibre direction surrounded by layers of ECM, whose orientation is locally characterised through the surface unit normal vector . While we note that at the fibre boundary, these layers are perfectly tangential to the individual fibres such that ⋅ = 0, this might not be the case with increasing radial distance from them due to the competing influence of neighbouring fibres. Therefore, we introduce a vector ′ that lies in the plane of the ECM layer such that ′ ⋅ = 0, and coincides with at the muscle fibre and ECM interface. This vector is computed from yet another steady-state diffusion analysis with boundary conditions defined as follows: An insulated interface condition between muscle fibres and ECM was again introduced choosing material properties as MF ∕ ECM ≪ 1 and a non-zero difference 1 − 2 between the upper and lower surfaces of the cube is assumed so that The vector ′ is then defined from the computed flux as For later use, a second in-plane vector = × ′ is introduced (Fig. 5) such that the triad of vectors { , ′ , } forms a local orthonormal basis in the ECM layer.

Reducing boundary effects
The solutions to the aforementioned BVPs are notably influenced by the boundary of the DOI which introduces unintended boundary artefacts. For instance, portions of muscle fibres along  are lost when the muscle fibre cube is sectioned out of a larger assembly (Fig. 4). To address this, a boundary extension to the DOI  is introduced (Fig. 6a) by selecting a larger hollow cuboid region enclosing the DOI (cf. Section 2.4). Correspondingly, surfaces , ECM , MF with = t, b, l and IF are defined (cf. Eq. (6) and Fig. 6e). We discuss this in 2D in a parametric study presented in Appendix B, wherein the effect of the size of boundary extension on iso-contour determination inside the DOI is explored.

Preliminaries
In this study, both muscle fibres and ECM are modelled as hyperelastic solid bodies wherein each material point of the DOI with position vector in the reference state is mapped onto the position vector = ( , ) in the current state at time t, with the deformation gradient ( , ) = Grad ( , ). The strain-energy density function of each component depends on through the right Cauchy-Green tensor = T . The presence of muscle fibres is taken into account by the use of an anisotropic invariant (Spencer, 1984), while the structural disposition of collagen fibres about mean preferential directions in the ECM layers is considered through the generalised structural tensor approach Holzapfel et al., 2015).

Muscle fibres
The behaviour of muscle fibres is modelled through a modified form of the transversely isotropic exponential strain-energy function proposed in Holzapfel et al. (2000) as where 0 , 1 , 1 and 2 are the material parameters (Holzapfel et al., 2000), and ⟨⋅⟩ denote Macaulay brackets. Here, C = tr denotes the first principal invariant of , = √ det represents the volume change, and M = ⋅ is the squared stretch in the local muscle fibre direction in the reference state. In contrast to (Holzapfel et al., 2000), a compressible form of the neo-Hookean type strainenergy function is chosen to model the isotropic behaviour of muscle fibres in Eq. (13) (cf. (Stracuzzi et al., 2018)). Muscle fibres contain a substantial amount of water and hence are commonly modelled as (nearly) incompressible Rehorn et al., 2014). Here, the parameter 0 > 0 controls the resistance to volume changes and thus allows to relax this assumption.

ECM layers
The architecture of the ECM as described in Section 3.1 features similarity with other quasi-planar tissues, such as skin (Limbert, 2017;Annaidh et al., 2012;Meyer et al., 1982), cornea (Pandolfi and Vasta, 2012;Pandolfi and Holzapfel, 2008), or arterial walls (Qi et al., 2015;Holzapfel et al., 2015), and therefore similar strategies for describing this structure in terms of models stand to reason. Particularly, the use of generalised structural tensors, resulting from the integral of rank-one structural tensors over bi-variate spatial orientation density functions, has been proposed to capture this type of anisotropy Pandolfi and Vasta, 2012;Annaidh et al., 2012) and we make use of a similar approach here.

Characterisation of collagen fibril architecture
In an ECM layer with a local surface normal , the two mean directions 1 and 2 about which collagen fibres are dispersed are given with respect to the local orthonormal basis { , ′ , } by where denotes the angle between ′ and the mean direction , = 1, 2 in the plane spanned by ′ and . The 3D dispersion about these two main directions is defined through an orientation density distribution function ( ; , ) centred along , and generalised structural tensors are then defined as Holzapfel et al., 2015) where 2 denotes unit sphere. The local average squared stretch of fibrils at a location within the network of dispersed fibres is given by The orientation distribution function can be explicitly expressed in terms of the polar angle = arccos( ⋅ ), ∈ [0, ] and the azimuth angle ∈ [0, 2 ] between and − cos( ), and is represented as a product of two distribution functions as  ( ; , ) =̃( , ) = ip ( ) op ( ).

Strain-energy function
To model the mechanical response of ECM, the generalised structural tensor approach is employed in conjunction with a hyperelastic variant of the soft tissue model proposed in (Rubin and Bodner, 2002). Thus, we obtain the strain-energy density function The terms 0 ( ) and S ( , 1 , 2 ) therein read 0 ( ) = 1 ( C − 3) + 1 2 0 ( −2 0 − 1) and R. Kuravi et al.

Material parameters
Material properties for muscle fibres are deduced from experimental data on quasi-static uniaxial tension and compression of single fibres . To this end, an analytical form of the material model described in Section 4.2 was used in conjunction with a custom Matlab script to identify the material parameter set  = { 0 , 1 , 1 , 2 } (Table 4). In Fig. 7a, the parametrised analytical model is compared with the experimental data.
To validate the calibrated model, particularly with regard to the compressive properties, transverse compression as described in Böl et al. (2019) was simulated based on a FE model generated from a crosssectional image of a single muscle fibre presented in Böl et al. (2019), and discretised with 5570 C3D4 type elements (Fig. 7c). Transverse compression between two plates was simulated with Abaqus/Standard software (Abaqus 6.14-1, Dassault Systèmes ® Simulia Corp., Providence, RI, USA) using rigid surfaces in contact with the fibre model, assuming hard friction-less contact. The predicted transverse compression of the single muscle fibre in Fig. 7d shows sound agreement with the experiment .
To the best of the authors' knowledge, there is currently no experimental data available on the large strain anisotropic behaviour of isolated muscle ECM, besides a single set obtained from uniaxial tension of de-cellularised muscle samples (Gillies et al., 2010). We are currently establishing a corresponding protocol, and first results were indicative for the range of the expected response. The material parameters of the ECM model (Section 4.3.2, Table 4) were thus identified such that a reasonable agreement with preliminary tensile tests was achieved. We note however, that the aim of this parametrisation is merely a match of the orders of magnitude with experimental data at this point. The Table 4 Material parameters identified for the passive response of muscle fibres, and the ECM.

Material
Parameter Value Muscle fibres

FE implementation and work-flow
Thermal and mechanical analyses were performed through a FE simulation of the corresponding BVPs in Abaqus/Standard software. The three BVPs of steady-state diffusion described in Section 3 were implemented and solved as thermal analyses on the extended domain  ∪ encompassing both the cubic sample () and a boundary extension () with their interfaces connected through a tie constraint (see also Appendix B). The element-level heat flux vectors obtained from each case were then used to define the local directions of anisotropy { } in the muscle fibres (Eq. (8)) as well as the local iso-surface normals { } (Eq. (10)) and the vectors { ′ } in the ECM (Eq. (12)) for every element .

R. Kuravi et al.
For mechanical analyses, material models for muscle fibres and the ECM described in Section 4 were provided as user-defined material subroutines (Abaqus/Standard 6.14-1, 2014b) in which the components of the vectors ′ , , and , obtained from aforementioned thermal analyses, were introduced as state variables in the undeformed configuration.
Finally, various mechanical BVPs were analysed on the DOI (Section 5.3), and the homogenised stress response was computed (Section 5.2). The entire work flow is summarised in Fig. 8. Table C.5 summarises the details of the FE discretisation for virtual muscle samples generated from one, two, and five images.

Determination of homogenised response
The effective response of the DOI was estimated by subjecting it to affine displacement boundary conditions and using averaging theorems in finite kinematics (Hill, 1972).
To this end, the volume averaged Cauchy stress tensor ⟨ ⟩, henceforth referred to as effective Cauchy stress, under quasi-static loading conditions with negligible contribution from body forces is obtained as (Costanzo et al., 2005) where and , respectively, represent volume and position vectors in the current state and̂represents the surface traction vector. In the FE formulation, the affine boundary conditions are realised by displacing all boundary nodes as ( ) = ( 0 − ) + on , where and , respectively, represent the displacement and the reference position of node on the boundary, while 0 , , and , respectively, denote the 'macroscopic' homogeneous deformation gradient, the identity tensor, and a constant rigid translation vector. Finally (Geers et al., 2010), represents the discrete form of ⟨ ⟩, where defines the total force at each node on the boundary. The volume averaged first Piola-Kirchhoff stress, henceforth referred to as effective nominal stress, is readily obtained as (Costanzo et al., 2005) We note that is set to zero, without loss of generality, by conveniently aligning the centroid of the DOI with the origin.

Load cases for mechanical analyses
The mechanical responses of the muscle samples to uniaxial extension with lateral contraction (UAE), uniaxial compression with lateral expansion (UAC), and simple shear loading were studied along different directions. Moreover, to illustrate the potential effect of 3D structure on active tissue properties, a ''pseudo-activation'' load case was considered by pre-stressing the muscle fibres and solving the global boundary value problem. To define the load cases, a set of orthonormal basis vectors { 1 , 2 , 3 } was introduced such that 3 is aligned along the gross muscle fibre direction, i.e. the stacking direction, along which the samples were reconstructed from histology.
We first consider UAE and UAC (i) along, (ii) transverse, and (iii) oblique (45 • ) to the gross muscle fibre direction, specified through deformation gradient  (Figs. 9j, k). The deformation gradient considered for simple shear reads where , ∈ {1, 2, 3}, { 1 , 2 } are normal to the lateral faces and specifies the amount of shear varied from −0.3 ≤ ≤ 0.3. Since symmetry of the shear response about = 0 cannot be assumed owing to the spatial heterogeneity of the muscle tissue, both positive and negative shear are considered for each of the 6 load cases (Figs. 9f-k). For each case of R. Kuravi et al. UAE and UAC, one of the principal stretches was predefined through a tensile stretch and the resulting two lateral stretches were numerically computed. Note that due to physiological asymmetry of the sample, two transverse directions and two 45 • directions were considered for the load cases (ii) and (iii). It is noted that due to heterogeneity and asymmetric material distribution, UAE and UAC load cases do not lead to the states of uniaxial tension and compression, and hence non-zero global shear stresses may occur.
The pseudo-activation load case is considered representative of a typical isometric tension experiment on a stretched muscle. It comprises of passive UAE up to = 0 = 1.3 along 3 followed by an isometric activation type increase of the muscle fibre stress. The latter is achieved by adding an active stress act to the axial fibre stress component (Röhrle et al., 2008), so that the Cauchy stress implemented in the user subroutine reads where act = with a predefined active peak stress 0 = 100 kPa (Fig. 9l, m), that represents a physiologically reasonable value (Edman, 1999;Teran et al., 2003).

Results
In the following, the results of the numerical simulations performed on a muscle sample reconstructed from 5 images (Sections 6.1 and 6.2) and a comparison of the response of samples developed from 1, 2 and 5 images, respectively, subject to the different load cases (Sections 6.3 and 6.4) are presented.

Computational estimates of anisotropy
Numerically estimated local directions of anisotropy in  MF and  ECM , i.e., the local muscle fibre vectors and the iso-surfaces are shown in Figs. 10a-d. As expected, the former (Figs. 10a, b) are fairly parallel to the gross direction of the muscle fibre sections, while the latter (Figs. 10c, d) are tangential and closely wrap around their surfaces.

Mechanical response of virtual muscle sample
The effective first Piola-Kirchhoff stress response and the corresponding volume changes of the 5 image sample in UAE/UAC (Figs. 9ae) is depicted in Figs. 11a,c. With the current parametrisation of the model, loading along the stack axis 3 yields the stiffest response, while the 45 • and cross-fibre responses are similar, although the latter are slightly softer (Fig. 11a). When comparing these predictions from one specific DOI to the available experimental data on porcine muscle tissues (Van Loocke et al., 2006;Böl et al., 2012), it is observed that the predictions lie within the typical order of magnitude of the expected stresses observed. However, the typical order of stiffness between tests in different directions is not captured (Fig. 11b). With the current set of parameters, the model also predicts notable volume changes in the UAE and UAC tests (Fig. 11d), which are most pronounced for tensile loading in cross-fibre direction, where they exceed 5% for a longitudinal stretch of = 1.3. Due to the spatial heterogeneity of the DOI and its missing symmetry with respect to the principal axes of deformation the homogeneous applied boundary conditions lead to non-zero (effective) shear stresses as expected. Therefore, Fig. 11c depicts plots of the maximal Cauchy shear stress for each of the UAE/UAC load cases considered. For instance, in along-fibre UAE/UAC, the largest shear stress component was 32 and for 45 • it was 31 .

Comparison of the mechanical response of 1, 2, and 5 image samples
The numerically simulated effective nominal stress responses of the 1, 2, and 5 image samples under UAE/UAC, respectively, are depicted in Figs. 12a-c. The results indicate only moderate differences between the extruded (1 image) and the other samples. For instance, comparing the 1 and 5 image models reveals that the latter predicts about 9%-12% lower stresses for a stretch of 1.3 along the fibre and 45 • directions, while the predicted stress for cross-fibre extension is about 23% higher. Likewise, the ascending/descending order of the loading directions in terms of the stress remains unchanged across the samples: While the along-fibre direction is the strongest in UAE/UAC, the 45 • and the cross-fibre direction exchange the next place from UAE to UAC. A closer inspection of the UAE/UAC reveals the dissimilarity between the two cross-fibre directions, and two 45 • directions, thus pointing to a lack of transverse isotropy. These differences are more pronounced in case of the extruded sample.
A larger effect of the 3D details encompassed in the 2 and 5 image samples can be observed in the response to simple shear, shown as effective (maximum) Cauchy stress in Figs. 12d-i. From Figs. 12d-f, it is observed that, in contrast to extruded (1 image) sample, 2 and 5 image samples exhibit asymmetric shear stress response to positive and negative shear in 31 and 32 modes (Figs. 9f, h). This becomes particularly evident in the resulting Cauchy stress normal to the shear plane as shown in Figs. 12g-i. This effect is attributed to the intrinsic tissue asymmetry about the plane normal to 3 in 2 and 5 image samples (Fig. 4j). Moreover, from Figs. 12d-f unequal shear stress response is observed between modes (i) 31 and 32 , (ii) 13 and 23 , and between (iii) 12 and 21 in all samples. These differences reflect the lack of material symmetries that would allow for rotations by ∕2 about gross muscle fibre direction 3 , such as transverse isotropy (see also Section 7.3). Interestingly, the difference between the shear responses of case (iii) becomes smaller with increasing 3D details, whereas cases (i) and (ii) show an increasing tendency from (d) to (f).

Pseudo-activation
The numerically simulated effective nominal stress response of 1, 2, and 5 image samples for the pseudo-activation load case described in Section 5.3 is depicted in Fig. 13a. The corresponding lateral stretches in the two transverse directions ( 1 and 2 ) normal to the faces of the muscle sample are shown in Figs. 13b, c, and Fig. 13d shows the corresponding changes in the volume. A movie exemplifying the response of the 1 and 5 image samples is provided as supplementary video material. The large differences between the 1, 2, and 5 image samples becomes evident in the kinematic response of the samples after activation of the muscle fibres at = 0 : While the extruded (1 image) sample maintains its shape when the straight muscle fibres stiffen, the 2 and 5 image samples exhibit unequal lateral dimensional changes resulting from the misalignment between local muscle fibre and stacking directions. Along with these changes, the 2 and 5 image models also predict a slight volume loss upon isometric activation with the given set of parameters.

Discussion
The goal of this paper was to evaluate the role of the 3D microstructure, in particular its variation along the gross direction of the fascicles, on the mechanical behaviour of skeletal muscle tissue. To this end, a method was developed to (i) segment stacks of histological sections obtained from muscle tissues into muscle fibre and ECM fractions, and (ii) to reconstruct 3D models of the tissue suitable for FE analyses.

Reconstruction of the internal 3D structure
A thorough investigation of the histological sections of a stack reveals notable variations in both the shape and the position of the cross-sections of a muscle fibre through the stack (cf. Farrell and Fedde, 1969). Therefore, the 3D models generated from more than 2 images, in this paper exemplified by the 5 image sample, are characterised by undulant fibres of non-homogeneous cross-section including those fibres that start or terminate within the stack. To ensure that these observations do not result from preparation artefacts caused by the cutting and staining protocols or errors in image registration and segmentation, we hypothesised that the curvy shape should also be evident in longitudinal tissue sections. To this end, a few additional longitudinal histological sections were prepared according to the established protocol. The examples shown in Figs. 14a-c reveal three major observations: (i) The longitudinal axes of muscle fibres generally exhibit a small variation along the fibre length (e.g. locations A, B and C in Fig. 14a). (ii) The cross-section area of muscle fibres varies along their length (e.g. locations A-C in Fig. 14a, and I-K in Fig. 14b). (iii) Several muscle fibres appear to either exit or terminate at the given histological section (e.g. arrows in Fig. 14a and c). These observations, though sample specific, confirm the non-uniformity along the muscle fibres (Farrell and Fedde, 1969), and underline that a sample extruded from a single histological section can generally not adequately represent the tissue scale micro-structure of skeletal muscle.
From the generation of 3D virtual samples ( Fig. 4d-g) it is observed that with increasing number of images chosen, the 3D structure of the muscle fibres gains in complexity including fibres that appear or vanish inside the muscle sample. We have attributed such fibres to intrafascicularly terminating fibres (cf. Sharafi and Blemker, 2011) and, noteworthy, for the muscle samples virtually generated, this concerns merely 4-5 of all accounted fibres (≈ 350).
Muscle fibre statistics deduced from the segmented image stack  reveal nearly 8% reduction in fibre continuity fraction (28 fibres) between images 1 and 13 , though only a smaller variation is observed in the total number of muscle fibres per image (Table 3). The closer analysis shows that this difference is largely concentrated along the boundary of the images where fibres enter or leave the images and appear as partial cross-sections. This could be attributed to the waviness of the fibres and their non-alignment with the stack axis.
In the present work, we used microtome-cut serial histological sections as a basis to generate morphologically realistic models of muscle tissue. Compared to non-destructive methods such as confocal or multiphoton microscopy (Nakamura et al., 2007;Campagnola et al., 2002), serial histology, even though it requires image registration (Section 2.2.1), guarantees similar image quality along the stack direction. In addition, it poses little restriction on the thickness of the samples to be reconstructed whereas in laser microscopy, the imaging depth depends on the optical density of the tissue, which is limited to about 100 μm (Plotnikov et al., 2006) for muscle, unless clearing agents are used (e.g., Plotnikov et al., 2006;Decroix et al., 2015;Williams et al., 2019). However, clearing or contrast-enhancing preparations, that also render muscle tissue amenable to micro computed tomography (Schaad et al., 2017) or serial block-phase scanning electron microscopy (Gillies et al., 2014), can cause tissue deformations such as shrinkage (Schaad et al., 2017) or affect the integrity of tissue components (e.g., collagen Plotnikov et al., 2006), thus manifesting as artefacts in the reconstructed geometry. We, therefore, consider serial section histology as a relatively simple and well-established method to reconstruct muscle samples up to the length scale of several millimetres. Notwithstanding, the methods presented here in Sections 3 and 4 could likewise be applied to tissue morphologies generated from alternative microscopy imaging techniques.

Computational estimation of anisotropy
The internal micro-structure of the individual components, i.e. muscle fibres and ECM, was below the length scale identifiable from the histological micrographs. However, this micro-structure, responsible for the anisotropic properties of these components, strongly affects the tissue scale response. To this end, computational estimates of component specific anisotropy were obtained by solving steady-state heat transfer problems, and the corresponding heat flux vectors were utilised to obtain principal directions of anisotropy.
Similar approaches were established, for example, to estimate the course of fascicles through a muscle (Choi and Blemker, 2013;Handsfield et al., 2017), or to identify the layered collagen architecture in arteries (Raina and Miehe, 2016). While our results lack a quantitative comparison with microscopy, they are in line with recent confocal microscopy studies on chicken muscle  revealing a layered ECM structure surrounding fibres and fascicles.

Effect of the internal 3D structure on the mechanical response
Until now, FE models of muscles have only been constructed from a single, either hypothetical (Spyrou et al., 2019) or histological (Marcucci et al., 2017;Blemker, 2010, 2011) section of muscle tissue, by extruding it perpendicular to the cross-sectional plane. This approach was incorporated as a special case in this study by considering a single image of the stack whose mechanical response is shown in the first column of Fig. 12. Comparing this with the response of models generated from 2 and 5 images reveals that though differences exist, their magnitude is highly dependent on the load case. While UAE and UAC load cases reveal moderate differences, the shear response, in particular the corresponding normal stresses perpendicular to the shearing direction, differ strongly. These results showcase the material asymmetry in the samples that is caused by the change in position and shape of the fibre cross-sections from one histological section to the next, which in turn manifests as the non-alignment angle of the  stack axis (here 3 ) with the muscle fibre directions in the DOI. This becomes evident from the mean fibre angles (Table 3) calculated from the muscle fibre centres with respect to the stack axis. For instance, in the 2 image sample, the stack axis deviates from the gross muscle fibre direction by ≈ 12.53 • on average. Unless the fibres lie in the shear plane, positive or negative shear will therefore set the muscle fibres initially under tension or compression, and thus lead to the observed asymmetry (Figs. 12e,f,h,and i). This non-alignment angle has major contribution from the sample preparation if the cutting, i.e. the cross-sectional plane is not perfectly perpendicular to the muscle fibre direction. However, this systematic error can hardly be mitigated since histological sections are obtained from much larger muscle samples and the muscle fibre directions vary locally. In this work, we did not rectify this misalignment by a corresponding rotation because, in fact, one would not know this angle if provided with a single histological section. This angle thus belongs to the systematic errors of using a 1 image model and should hence be included for a fair comparison. Nevertheless, the main effect of this misalignment can be estimated from the difference between the 1 image (extruded) and the 2 image model, which consists of straight but tilted fibres. Noteworthy, even the 1 image sample response reveals deviations from symmetry that allows rotations by 90 • about the stack axis, and in particular transverse isotropy, in line with the results on fascicle scale models in Sharafi and Blemker (2010). Therefore, at the scale investigated here, the material symmetry of skeletal muscle tissue seems more complex than transverse isotropy, and may lead to local shear deformations within the tissue. These observations are of practical importance, particularly towards the development of microstructure driven continuum constitutive models for muscle tissue (e.g., Gindre et al., 2013;Bleiler et al., 2019), because local shear is practically unavoidable and lateral force transmission in muscle through the ECM shear has been hypothesised to play an important role Blemker, 2011, 2010).
Another noteworthy effect of including 3D structure is a slight reduction in the well-known tension-compression asymmetry in the uniaxial mechanical response of skeletal muscle tissue , observable in the fibre direction response red curves). This can be explained by the existing 'undulations' of muscle fibres and the ECM (along the gross muscle fibre direction) that attenuate the strong change of stiffness between tensile and compressive loading as compared to the extruded model, where fibres are straight so that the reference state marks the tension-compression transition.
The response of 1, 2, and 5 image samples to the pseudo-activation load case shown in Figs. 13a-d exhibits the most revealing influence of the internal 3D structure. Here, substantial stiffness is added along the local muscle fibre directions upon activation. Hence, to minimise the energy, the 2 and 5 image samples change their configuration (Supplementary Video) so that the local fibre vectors deform, and the internal loads are redistributed among the muscle tissue components. This is in contrast to the extruded model, wherein no further dimensional changes occur since the straight fibres are already in the energetically optimal state with respect to the given, simple isometric activation load case.
In addition to studying the effect of 3D microstructure on the mechanical behaviour of muscle tissues, the here presented 3D models could be used to perform virtual experiments for various other computational studies on skeletal muscle tissue. Potential research questions concern, e.g., interactions between individual components, for instance, in relation with tension-compression asymmetry , the effect of variation of properties arising due to muscle activation or muscle pathology, e.g., DMD (Virgilio et al., 2015) or the development of microstructure driven continuum constitutive models for muscle tissue (e.g., Gindre et al., 2013;Bleiler et al., 2019). To this end, selected models are provided as additional supplementary material (Kuravi et al., 2020).

Predicted stress response
The effective tensile stress responses of the 5 image sample (Fig. 11a) were obtained after parametrising the constitutive models only from single muscle fibre data and preliminary experimental results on ECM in uniaxial tension tests. The generally sound agreement of the predicted tissue response with existing uniaxial compression data on porcine muscle tissue (Van Loocke et al., 2006;Böl et al., 2012) in terms of order of magnitude suggests that the model is able to capture key mechanisms responsible for the tissue scale behaviour. Yet the model predicts that in UAC, along-fibre direction provides the overall stiffest response, followed by the cross-fibre and 45 • directions (Figs. 11a,b). This is different from the order observed in the experiments (Van Loocke et al., 2006;Böl et al., 2012), which reveal the highest stresses in the cross-fibre direction, followed by the 45 • and fibre directions (Van Loocke et al., 2006) or vice-versa (Böl et al., 2012). The potential causes for these deviations are manifold, including but not limited to (i) the experimental protocols, (ii) muscle type and its age, and (iii) it should be noted that the DOI used for the simulation represents one specific tissue sample, whereas the experimental results are the mean values obtained from a large set of tests conducted on much larger samples. Notwithstanding, the altered order observed in the simulated response curves in terms of stiffness, in contrast to the experiments, could be attributed to the material parameters of the anisotropic nonlinear ECM material model (Section 4.3) deduced from yet incomplete experimental data, see also Section 7.5.

Limitations
The first step of image registration consists of an affine transformation of the stack  . This operation is based on two key assumptions: (i) the first image is undeformed and can be considered as a reference against which the subsequent images need to be adjusted by alignment of the pin-hole markers, and (ii) that a single global affine transformation is sufficient to capture all deformations induced by the cutting process. Both assumptions could be resolved by using more pin hole markers and calibrating their initial foot print.
Following the affine registration, a clustering-driven segmentation of the images is performed. This operation is influenced by the number of clusters used. Following preliminary trials, a choice of 6 clusters to isolate muscle fibres from the image stack and a choice of 4 clusters for low level clustering were observed to be adequate (Table 1). This choice, however, depends on the effectiveness of the staining protocol, the image quality and the operator. Despite structural differences between endomysium and perimysium, and the presence of other components such as fat, only muscle fibres and collagenous ECM were chosen as primary components. Nevertheless, with aptly chosen staining protocols and clustering parameters one could further include other components.
In view of the scope of the present paper, high confidence in the matching of muscle fibre cross-sections belonging to the same fibre in different images of the stack is vital. Their identification through an operator was therefore preferred over an automatic procedure. In this way, problems of non-uniqueness associated with an automated procedure can be avoided. In fact, by an operator choice we could guarantee that fibres would neither intersect nor bifurcate in the virtual model, pointing at a unique mapping. To facilitate this, the use of fascicle and fibre details from additional slices, that were not included in the reconstruction process, served well to track the fibres through the relevant images.
The availability of experimental data for individual muscle components is extremely limited, and even less is known about the interactions between them. Due to this limitation, we assumed that the displacements of ECM and muscle fibres are fully coupled at their interface, although the simulation framework would generally allow for various other contact interactions. While the dedicated connections between the collagen fibrils in the ECM and muscle fibres reported from SEM images (Purslow, 2008;Passerieux et al., 2007) point at potentially intricate interaction mechanisms, the representation of ECM as a continuous material in our simulations does not allow to account for connections formed by single collagen fibrils. The inadequate parametrisation and validation of the ECM material model (Section 4.3) is currently one of the main limitations regarding its predictive abilities. Due to the complex material symmetry, several independent experiments would be required for calibration. Such data is, to the best of the authors knowledge, not available to date and will be the subject of future work of our groups.
The analysis presented in this paper is based on models generated from a single histological stack of images reconstructed with great level of detail, although the methods were tested for other samples as well. For a statistical analysis of these results, a large number of samples would be needed. The large number, in turn, would call for a precise automated image registration and segmentation procedure for efficient implementation, which is beyond the scope of this paper.
Finally, it should be noted that the DOI, chosen here as 300 μm × 300 μm × 300 μm might be too small to be considered as a representative volume element of the tissue. This size was selected as a compromise between precise representation of tissue-scale properties and computational efficiency. Since continuum models of skeletal muscle typically homogenise on the tissue and not fascicle scale, the DOI was chosen to be large enough to contain not only endomysium but also substantial amounts of perimysium.

Summary and conclusions
In this work, an approach towards developing 3D virtual muscle tissue samples of desired size and level of detail from stacks of histological sections was presented. The tissue components considered for modelling were limited to muscle fibres and the extracellular matrix, but the latter could be more refined to incorporate other muscle tissue components, such as fat, by adapting the staining and segmentation protocols. Through a combination of image clustering and Voronoi tessellation methods, the proposed procedure is able to yield realistic volume fractions and spatial distributions of muscle fibres and the extracellular matrix. The thus established 3D computer models of muscle tissue were discretised with finite elements, and equipped with non-linear anisotropic hyperelastic material properties. A steady-state diffusion based method was developed to estimate the local directions  Fig. B.16. Effect of the ratio of edge-lengths of the  ∪ and , i.e., on average angular deviation ( ), for a numerical tolerance of = 10 −6 . of anisotropy that were below the resolution of the histological micrographs. These computer models were used to quantify the importance of capturing the change of micro-structure along the gross direction of fascicles, i.e. the slightly wavy and inclined course of the muscle fibres through the tissue, and the non-uniform spatial disposition of the extracellular matrix. Finite element simulations performed with models generated from 1, 2, and 5 images clearly show that at the length scale of a few fascicles, here 300 μm, muscle tissues are not transversely isotropic, owing to asymmetric material distribution not only about the gross fascicle direction but also with respect to the planes perpendicular to it. More precisely, the micro-structure was shown to cause high heterogeneity in deformations and stress within the tissue even under homogeneous boundary conditions. The responses to shear loads and muscle activation wherein the muscle fibres drastically increase their stiffness, are particularly affected by this, while relatively less influence was observed for tensile and compressive load cases. This indicates that finite element models generated by extruding a single histological section cannot adequately capture the muscle response to generic loads, while they may still apply to special cases. In conclusion, our study strongly suggests to incorporate the internal 3D micro-structure in multi-scale finite element models, and to revise the assumption of transverse isotropy in the corresponding continuum mechanical constitutive models of skeletal muscle tissue.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Figs. B.15b and B.15c, respectively, depict the isocontours inside  ECM for = 1.0 and = 1.8. In both cases, it is observed that away from the boundary, the isocontours bear resemblance to the layers of ECM around the muscle fibres . However, in the former case, the boundary effect (Fig. B.15c) occurs which is quantified through . From Fig. B.16, it is observed that , in general, decreases with . It is observed that for = 1.667 (i.e., = 500 μm), the angular deviation yields = 1.74 • , i.e., a reasonable compromise between the boundary effect and the computational time. This corresponds to adding about one typical fibre diameter on each side of .

Appendix C. Mesh dependency
The choice of the mesh size for the FE analyses is validated by refining the mesh used in FE analyses with approximately four times more elements (Table C.5). The along-fibre UAE case is considered as a representative case and the outcomes from these two meshes are compared in Figs. C.17a-d. Fig. C.17a depicts nominal along-fibre stress, and Figs. C.17b-d, respectively, depict shear stresses 12 , 13 and 23 for both the mesh sizes. The negligible difference between the outcomes of these two dicretisations suggests that the mesh size chosen is small enough for the FE analyses.

Appendix D. Supplementary data
Supplementary material related to this article can be found online at https://doi.org/10.1016/j.jmbbm.2020.104109.