Mechanical Biomedical Materials Predicting muscle tissue response from calibrated component models and histology-based finite element models

rigorously predictive and thus unforgiving strategy suggest that the prediction of the tissue response from the individual characteristics of muscle cells and decellularised tissue is only possible within clear limits. While orders of magnitude are well matched, and the qualitative behaviour in a wide range of load cases is largely captured, the existing deviations point at potentially missing components of the model and highlight the incomplete experimental information in bottom-up multiscale approaches to model skeletal muscle tissue.


Introduction
Mathematical and computational modelling of skeletal muscle tissues poses a special challenge due to the composition, complex spatial distribution and orientation of its constituents, and their respective internal hierarchical structure, which altogether contribute to the overall motion, force generation and gait stabilising characteristics of skeletal muscle. The tissue is predominantly comprised of long, multi-nucleated excitable muscle fibres grouped into fascicles to finally form muscles (Barrett et al., 2016;Lieber, 2002). This hierarchy is established by collagenous layers referred to as endomysium, perimysium and epimysium that wrap fibres, fascicles and muscles, respectively (Barrett et al., 2016;Lieber, 2002). These layers form a substantial part of the extracellular matrix (ECM) and constitute 1%-10% (Kjaer, 2004) of the muscle dry weight. The layers vary in their composition, structure and also in their respective mass fractions. For instance, endomysium exhibits a thin continuous random network-like structure shared by adjacent offer ease of numerical implementation and provide an 'averaged' macroscopic response, they are unable to relate this response to internal changes of microstructure, to the distribution of stresses and strains among components, or to their interactions. In this regard, the use of multi-scale models encompassing the microstructure, with distinct continuum representations for the individual constituents, i.e. muscle fibres and ECM, have recently gained importance to explore the effect of microscopic properties on the macroscopic behaviour (Sharafi and Blemker, 2010;Bleiler et al., 2019;Gindre et al., 2013), and to address medical questions in relation to muscle pathologies such as Duchenne muscular dystrophy (DMD) (Briguet et al., 2004;Emery, 2002), spasticity (Lieber et al., 2004) and deep pressure ulcers (Stekelenburg et al., 2006). While tension-compression asymmetry was explored in Gindre et al. (2013) using a structural model, more involved finite element (FE) based micromechanical models were employed to explore the microstructure dependent longitudinal shear modulus (Sharafi and Blemker, 2010), the effect of DMD-related tissue alterations on the damage susceptibility (Virgilio et al., 2015), and the component specific contributions to passive muscle behaviour (Marcucci et al., 2019). Multi-scale considerations were also utilised in conjunction with homogenisation frameworks to develop continuum constitutive models (Spyrou et al., 2017(Spyrou et al., , 2019Bleiler et al., 2019). However, these microstructure models are based on either information from single histological sections of muscle tissue cut perpendicular to the muscle fibres (Marcucci et al., 2017(Marcucci et al., , 2019Blemker, 2010, 2011), on artificially generated cross-sections (Virgilio et al., 2015;Spyrou et al., 2017Spyrou et al., , 2019 or on idealised tissue geometries in combination with the assumption of affinity (Gindre et al., 2013;Bleiler et al., 2019). Hence, these current approaches do not incorporate microstructural changes along the fibre axis, which break symmetry and add additional model complexity. In a recent study we have reported on the effect of including this complexity by comparing the responses of an extruded micro-scale FE model and full 3D models for various load cases .
While there are comprehensive experimental data sets on tissue scale passive behaviour of skeletal muscle, e.g. under uniaxial tension (e.g. Calvo et al., 2010;Morrow et al., 2010;Mohammadkhah et al., 2016;Takaza et al., 2013aTakaza et al., , 2014Nie et al., 2011;Gras et al., 2012Gras et al., , 2013, unconfined and semi-confined compression (Van Loocke et al., 2006;Böl et al., 2012;Van Loocke et al., 2009;Chawla et al., 2009;Morrow et al., 2010;Takaza et al., 2014Takaza et al., , 2013bVan Loocke et al., 2008;Böl et al., 2014Böl et al., , 2016, shear (Morrow et al., 2010;van Turnhout et al., 2005) as well as on single muscle fibres (Rehorn et al., 2014;Böl et al., 2019;Toursel et al., 2002;, the mechanical characteristics of ECM have remained relatively unexplored (Lewis and Purslow, 1989;Lewis et al., 1991;Purslow and Trotter, 1994), and even more, direct experimental information on the mechanical inter-component interactions is missing. Towards shedding light on the anisotropic mechanical behaviour of the isolated collagenous ECM, we have recently completed a comprehensive study on decellularised muscle samples (Kohn et al., 2021). In the present contribution these novel data on ECM material characteristics are supplemented by equibiaxial tests and incorporated in detailed micro-structural FE models  to investigate the anisotropic mechanical behaviour of skeletal muscle tissues. To this end, ECM is associated with a constitutive material model inspired from its physiology, and simulations on virtual, decellularised muscle samples are performed in order to estimate the related material parameters through inverse FE analysis. The calibrated ECM FE model is combined with the muscle fibres, and the response of the composite muscle samples under various loading conditions is simulated. These predictions are compared with existing experimental data and a new set of equibiaxial extension data on fresh porcine muscle. Finally, the model is used to study effects of sample size and the characteristics of the fibre-ECM interface.

. Ethical approval
The study was exempted from ethical committee review according to national regulations (German Animal Welfare Act), as porcine skeletal muscles of healthy, female domestic pigs were obtained from a slaughterhouse immediately after animal sacrifice.

Sample preparation
Hind legs ( = 2) of female domestic pigs (Sus scrofa domestica, age: ∼6 month, mass: ∼90 to 100 kg) were obtained from a slaughterhouse immediately after animal sacrifice. After the legs were transported to the laboratory in a cool box at 4 • C, biceps femoris muscles were excised and rectangular samples with characteristic edge lengths of approximately 60 mm and a thickness of 6 mm were cut. Care was taken to ensure that two edges of a sample were aligned parallel to the fibre orientation ( Fig. 1). Further, alginate was used to stabilise the tissue during cutting with a utility cutter (Böl et al., 2014). To prevent rigor mortis during the preparation and storage period, the fresh muscles and samples were wrapped in cloths, soaked with Dulbecco's phosphate-buffered saline (DPBS) at 4 • C in a climatic chamber.
Before performing equibiaxial extension tests ( = 14) the differently processed tissue samples were divided into two groups (G 1 , G 2 ) and the specimens were trimmed to edge lengths of approximately 40 mm (Fig. 1, Table 1). The actual edge lengths ( 1 , 2 ) in both directions were determined from digital images as mean values of two opposite sample edges and the thickness was defined as the mean value of the thicknesses of all 4 edges. G 1 comprises the fresh, untreated tissue samples as excised from the muscle, whereas G 2 includes NaOHtreated samples with muscle cells removed according to a recently presented protocol (Kohn et al., 2021). In brief, the fresh samples were placed in 3% NaOH solution for 18 h. Thereafter, the treated tissues were removed from the NaOH solution and placed for 45 min in DPBS solution and then again in fresh DPBS solution for another 15 min, for more details see Kohn et al. (2021).

Biaxial testing equipment
Mechanical experiments were performed on a planar biaxial testing machine (Zwick GmbH & Co. KG, Germany) with four independently controllable linear actuators (Fig. 2a). Forces were measured by two load cells along each axis, with maximum load capacity of 100 N and  Experimental setup: (a) View of the tissue specimen mounted in the testing machine and (b) idealised illustrations of the marker (filled cycles) and hook (filled squares) positions. Note, the dimensions are given in millimetres, and the directions 1, 2, and 3 of the coordinate system correspond to the transversal, longitudinal, and -direction, respectively (cf. Fig. 1). 1 mN resolution, coupled to mounting devices that each consist of a U-shaped aluminium profile and a bolt holding a set of five hooks and cords (Fig. 2a). Additionally, a video extensometer (Video Extensometer ME46, Messphysik Materials Testing GmbH, Austria) with a minimum resolution of 0.4 μm was mounted above the specimens for tracking auxiliary markers on the upper side of the specimen (Fig. 2b), and whose output is used to control the movement of the actuators. A monochrome CCD camera detects the light-to-dark transitions on the surface of the samples and produces digital pixel images in 256 shades of grey. The software used for the measurement system (videoXtens, Zwick GmbH & Co. KG, Germany) processes the real time data recording and enables the tracking of the markers. Load cell values, actuator positions, and marker displacements were measured at 20 Hz and used for data analysis (Section 2.1.5).

Sample processing and testing protocol
Each dissected square sample was fixed with 20 hooks (five per sample side) regularly distributed along the sample edges, with a 7.0 mm hook-to-hook distance and 5.0 mm corner hook distance (Fig. 2b), see Morales-Orcajo et al. (2018) for details. Black circular markers (diameter 1 mm) distributed in a nine-dot square grid (3 mm distance between the markers) were glued on the upper side of the specimen to enable in-plane displacement measurements with the video extensometer (Fig. 2). The tracking area constitutes 4% of the area bound by the hooks (30 × 30 mm), which ensures a homogeneous stress distribution in the tracking area (Sun et al., 2005). Finally, the prepared specimen was positioned on the mounting devices of the biaxial machine, so that the muscle fibre direction was aligned with 2 ( Fig. 2a). To remove the weight effect and ensure the correct measurement of the initial markers distance in-plane, for all tests a preload of 10 mN was applied in both axes. The stretch-controlled experiments were realised with a strain rate of 0.5% s −1 until sample failure. Due to the short test duration of a maximum of one minute, it was not necessary to moisten the sample during the experimental realisation.

Data analysis
Under the assumption that the deformation is homogeneous in the central region of the specimen, the two stretch ratios along and perpendicular to the visible main course of the muscle fibres were determined from the tracked displacements (e.g. Morales-Orcajo et al., 2018), and controlled to be equal in order to generate at state of equibiaxial extension, characterised through the stretch ratio . They represent the components 11 and 22 of the deformation gradient (cf. Section 2.3) with respect to the coordinate system given in Fig. 2.
R. Kuravi et al. Fig. 3. Summary of all the steps involved in the development of 3D virtual muscles from histological images according to Kuravi et al. (2021). Source: Images partly adapted from Kuravi et al. (2021). They are not principal stretches in general, due to the presence of shear in the heterogeneous, anisotropic samples.
The two perpendicularly directed forces 1 and 2 were determined as the averages of the two opposite force sensors, respectively, and the nominal stresses were calculated as (1) using the sample-specific dimensions ( 1 , 2 , ) of each specimen (Table 1). Linear interpolation between the dense data points was used in order to compute mean and standard deviation at representative values of stretch, and to represent them as nominal stress vs. stretch, i.e. − curves.

Generation of 3D FE models from histological sections
3D virtual muscle samples are either imported from existing data sets (Kuravi et al., , 2020 or generated from histological images as described in Kuravi et al. (2021). Briefly, the following steps are performed to achieve these models: A set of pre-selected images of the histological image stack are aligned with each other through affine registration utilising auxiliary markers. A stack of square-shaped subimages is cropped, subjected to a set of image processing operations with Fiji (Münch, 2019;Schindelin et al., 2012), and Voronoi tessellation is performed. This yields a segmented image stack wherein individual muscle fibres and collagenous areas are identified. Boundary coordinates of individual muscle fibres are matched across all the segmented images and smoothly connected to form 3D virtual muscle fibres. The pore space is finally filled and attributed to the ECM, and the components are combined to a discretised muscle sample , which was meshed with tetrahedral elements (Appendix , Table A.3) such that they share same interface nodes, using FE software (Abaqus/CAE 6.14-1, Dassault Systèmes). A graphical summary of the procedure is provided in Fig. 3.
The local anisotropy of the model within the muscle fibre (MF) and ECM sections is determined by solving boundary value problems (BVP) of steady state diffusion . Briefly, the boundary  of the muscle tissue sample  is subdivided into top, bottom and lateral faces (  ,  and  ), which are in turn partitioned into fibre and ECM portions. In addition, all the internal interfaces between muscle fibres and ECM and their corresponding outward normal vector are identified as  IF and IF , respectively. Fibres and ECM are equipped with auxiliary coefficients of diffusion, and by specifying concentrations or fluxes at the boundaries and interfaces, various BVPs are analysed and the flux vectors associated with the local anisotropy are computed. To mitigate the boundary effect, BVPs were solved on extended domains ( ∪), obtained from 200 μm larger histological sections, that include boundary extensions of 100 μm (approximately one typical fibre diameter) width at the lateral faces and a tie constraint was enforced between  and (Fig. 3). The computational solution of the BVPs, implemented as heat transfer analyses in Abaqus/Standard (Abaqus 2018, Dassault Systèmes) provide for each element within the muscle fibre a unit vector that defines the local muscle fibre direction,  and within the ECM a triplet of vectors { ′ , , }. ′ is defined such that it is parallel to at the fibre-ECM interfaces  IF , specifies the local unit normal to the layered ECM structures, and = × ′ completes these vectors to an orthonormal basis, see Kuravi et al. (2021) for details.

Kinematics and preliminaries
In line with the standard notation in continuum mechanics, let the configuration of a body  in the reference and current (at time ) states be denoted by () and (), respectively. Each material point of  corresponds with the positions ∈ () and ∈ (), which are linked through the mapping = ( , ). The deformation gradient ( , ), its determinant ( , ), and the right Cauchy-Green tensor ( , ) are defined through where the dependence on place and time is understood. For an anisotropic hyperelastic material characterised by reinforcing families of fibres whose direction is specified by the unit vector , = 1, 2, … , , the strain-energy density function (SEDF) must be a scalar isotropic tensor function of and = ⊗ , = 1, 2, … , (Truesdell and Noll, 2004, Sec. 11). Invoking the representation theorems for isotropic scalar valued functions of symmetric tensors (Truesdell and Noll, 2004), the SEDF can be given in the form, see e.g. Schröder and Neff (2003), Itskov and Aksel (2004), Ehret and Itskov (2007) where 1 = tr , 2 = 3 tr −1 , 3 = det = 2 are the principal invariants of , and where typically the list of arguments of̃is restricted to a subset of these invariants. For materials with fibres dispersed symmetrically around the directions, generalised structure tensors can be defined (Advani and Tucker III, 1987;Freed et al., 2005;Gasser et al., 2006). For the special case where the mean directions, say , = 1, 2, … , , lie in one plane with unit normal (Fig. 4) and the dispersion can be characterised through two von-Mises type distributions that specify the in and out-of-plane distribution, the generalised structure tensors can be represented as (Holzapfel et al., 2015) = where ip ∈ [0, 1∕2] and op ∈ [0, 1∕2] represent in-plane and outof-plane dispersion parameters, respectively. The SEDF can then be represented as and corresponding invariants can be defined such as that replace the invariants ( ) 4 in the list of arguments of̃in Eq.
(3). The Cauchy stress is finally obtained as Based on these frameworks, we defined constitutive laws for muscle fibres and ECM in our previous work . The former was obtained by modification of the model in Holzapfel et al. (2000) with regard to compressibility as Kuravi et al. (2021) where M = tr( ), ⟨ ⟩ = (| | + )∕2, 1 and 1 are material parameters with unit of stress, and 0 and 2 are positive, dimensionless constants. The ECM model is discussed and amended in what follows.

Modification of the model for ECM
Drawing motivation from the experimentally observed muscle ECM microstructure, its composition and distribution in the muscle tissue (Oshima et al., 2007;Gillies and Lieber, 2011;Kjaer, 2004;Purslow, 2008Purslow, , 1989Light et al., 1985), and its resemblance with the network structures found in other tissues such as skin, cornea and arterial walls (e.g., Limbert, 2017;Pandolfi and Vasta, 2012;Holzapfel et al., 2015), the ECM was modelled based on the generalised structure tensor approach in Kuravi et al. (2021). In particular two local mean fibre directions were considered, symmetrically inclined at an angle to the mean muscle fibre direction ′ , and lying within the plane with unit normal , that defines the orientation of the collagenous layers (Fig. 4a). This concept was combined with the hyperelastic variant of the soft tissue model introduced in Rubin and Bodner (2002), and the corresponding SEDF was given by Kuravi et al. (2021) where in view of Eqs. (5) and (7) = 1 ( 1 − 3) + 1 2 0 ( −2 0 − 1) The now available experimental data on ECM (Kohn et al., 2021) pointed at a significant influence of shear across the layered ECM structures, which is in line with the hypothesis that the transmission of lateral forces across muscle fibres in muscle tissues is largely facilitated though shearing the ECM (Sharafi and Blemker, 2011;Huijing, 1999;Purslow, 2008). Therefore, the previous SEDF  is supplemented by a term that depends on a simple kinematic invariant related to the shear across ECM layers.

Numerical simulations 2.4.1. FE implementation and effective stress
For mechanical analyses, the material models of muscle fibres and the ECM (Section 2.3) were implemented through user-defined material subroutines (Abaqus/Standard 6.14-1, 2014), wherein the local material symmetry defined by the vectors { , ′ , } is introduced as state variables in the undeformed configuration.
The effective mechanical response of the muscle sample in a given BVP (Section 2.4.2) was estimated using averaging theorems in finite kinematics (e.g. Hill, 1972). To this end, affine displacement boundary conditions were invoked in discrete form as where and denote the displacement and the reference position of node on the boundary , respectively. 0 represents the 'macroscopic' deformation gradient, is the identity tensor of second order, and is a constant vector representing rigid translation, conveniently set to . The volume averaged Cauchy stress tensor ⟨ ⟩ and the corresponding effective nominal stress ⟨ ⟩ (volume averaged first Piola-Kirchhoff stress) are given as (Costanzo et al., 2005;Geers et al., 2010) where is the current volume of the modelled domain while and respectively are the total force and position of node on the boundary in the current state. The stress (19) results from volume averaging over the tissue (and not only ECM) domain. For the decellularised samples, this apparent stress is thus lower than the effective stress in the ECM, and clearly the factor in between is the volume fraction of ECM. Nevertheless, given the experimental difficulty to determine the volume and cross-section area of the decellularised muscle samples, this apparent stress is used for adjustment to the experimental data (Kohn et al., 2021).

Load cases
The load cases considered in this study comprise uniaxial extension with lateral contraction (UAE), uniaxial compression with lateral expansion (UAC), biaxial extension with lateral contraction (BAE), and semi-confined compression (SCC), as illustrated in Figs. 5a-l.
UAE/UAC of muscle samples and ECM are simulated along the longitudinal (0 • , Fig. 5a), transverse (90 • , Figs. 5b,c), and oblique (45 • , Figs. 5e,f) directions with respect to the gross muscle fibre direction, aligned with 3 . BAE is considered by stretching equally along this and one orthogonal direction (Fig. 9d). Finally, three modes of semi-confined (plane strain) compression are considered (Böl et al., 2014), where compression is imposed either along the muscle fibres (mode I) or along one of the cross-fibre directions (modes II, III) while one of the transverse dimensions was held constant (Figs. 5gl). Noteworthy, due to the lack of symmetry we considered two cases for each mode, fixing one of the two lateral faces, respectively. The tensile/compressive stretches were specified along the loading direction(s) while stretches along the unconstrained perpendicular directions were numerically computed to achieve traction free lateral faces. The deformation gradient is thus given as where {̂1,̂2,̂3} is an orthonormal basis and { 1 , 2 , 3 } are the corresponding principal stretches. The basis {̂1,̂2,̂3} is aligned with the global basis { 1 , 2 , 3 } except for oblique UAE, where it is rotated by 45 • about 1 or 2 , cf. Kuravi et al. (2021). For UAE/UAC one of the stretches is prescribed, while BAE of muscle and ECM samples (Fig. 8d) is specified by prescribing 1 = 3 . To simulate SCC of muscle samples (Figs. 5g-l) one of the stretches is kept at unity, one is prescribed through the experimental compressive stretch, and the non-constant transverse stretch is computed. Noteworthy, as a result of the histological preparation, there is a slight but practically unavoidable deviation between the gross muscle fibre direction and the stack axis 3 when generating the models, cf. . Based on the computed local fibre directions for each finite element, the mean orientation of all muscle fibres (̄M F ) with respect to the stack axis of histological images was evaluated. For the model considered here this angle was 17.06 • . To simulate loading along and across the mean muscle fibre direction, this misalignment was allayed to < 0.01 • by rigidly rotating the cubic muscle sample with an orthogonal transformation that maps̄M F onto 3 , a vector of the orthonormal basis { 1 , 2 , 3 }. Please note that this rotation is not shown in Fig. 5 for the sake of clarity.

Material parameters and inverse FEM
The parameters (Table 2) of the muscle fibre model (9) were adopted from Kuravi et al. (2021), where they had been obtained by fitting the model to experiments on individual muscle fibres (Böl et al., 2019). Inverse finite element (iFEM) calculations were performed to identify the set of material parameters of the ECM (Eqs. (16) and (17)). Briefly, a set of optimal material parameters that minimises the difference between experimental and simulated responses is deduced in an iterative fashion through an optimisation algorithm driven FE analysis, see e.g. Böl et al. (2013), Takaza et al. (2013b), Böl et al. (2012) and Ahn and Kim (2010).
Experimental data for the following three load cases was adopted from Kohn et al. (2021): UAE of the decellularised muscle tissue (ECM) along the longitudinal (0 • ), oblique (45 • ), and transverse (90 • ) directions with respect to the mean muscle fibre direction (Fig. 6). In the underlying experiments (Kohn et al., 2021), the muscle fibres were dissolved, thus leaving voids and causing the ECM structure to partially collapse. To avoid the computational effort with the consideration of self-contact between internal surfaces of the ECM, numerical simulations were performed on the whole muscle domain  but the fibre region ( MF ) was modelled as a weak, compressible neo-Hookean material by setting 1 = 0.00155 kPa, 0 = 1.0 and 1 = 0 in Eq. (9).
The optimisation algorithm to identify the parameter set ( opt ) was implemented through MATLAB-based (Version 9.3, R2017b, The MathWorks Inc., Natick, MA, USA) control utilising the Nelder-Mead simplex search algorithm (Lagarias et al., 1998) embedded in MATLAB as fminsearch function, minimising the objective function Here, exp and sim ( ), respectively, represent experimental and (effective) simulated nominal stresses for a given parameter set at the th increment of = 10 deformation increments. Noteworthy, in line with the experiments in Kohn et al. (2021) the nominal stresses represent the apparent stress that is obtained when relating the force acting on the decellularised samples to the original cross-section of the fresh tissue sample before the muscle fibres are removed. This avoids the complex determination of the void area. Fig. 7 summarises the optimisation algorithm.
In this way, the parameter set = { 0 , , 0 , 1 , ℎ 1 , ℎ 2 , 1 , 2 } was optimised whereas, to reduce the number of unknowns, the structural dispersion parameters ip = 0.15 and op = 0.49 (Eq. (5)) were preselected as a constitutive choice. Likewise, the parameter 2 controlling compressibility was set constant ( 2 = 10 5 ) since experimental information on volume changes was lacking. Noteworthy, this parameter affects the material compressibility, not the compressibility of the porous structure generated by removing the muscle fibres from the ECM. The   optimisation was terminated when the decrease in value of the objective function was marginal over several consecutive iterations, and the simulated curves showed sound agreement with the experimental ones upon visual inspection.

Sample types and component interactions
At first all simulations on muscle tissue and ECM were performed with a cubic sample of 300 μm edge length generated from 5 serial histological sections and developed in Kuravi et al. (2021). This model is available at Kuravi et al. (2020).
In Kohn et al. (2021) the effect of laterally compressing decellularised muscle samples before their response in UAE was studied. This procedure was repeated in the present study to perform BAE tests. The lateral compaction leads to flat and thin, membrane-like samples. To mimic this effect in the computational models, flat ECM models were generated by a purely geometrical operation: One of the transverse dimensions of a 300 μm cubic sample was scaled down by = 0.1, resulting in 90% 'compacted' samples ( Fig. 5a), and the directions of local anisotropy were calculated from the modified vector triad where ∈ { , ′ , } denotes the corresponding modified vector of . These samples were then subjected to UAE and BAE boundary conditions as described above and illustrated in Fig. 8. Note that the stress was still considered in terms of the apparent nominal stress, i.e. force per unit undeformed cross-sectional area of the muscle sample before decellularisation.
To assess the effect of (i) muscle sample size and (ii) muscle fibre-ECM interface properties, UAE/UAC and SCC load cases were simulated on thin samples based on a single histological section with a crosssection of 900 μm × 900 μm and a (physically irrelevant) thickness of 20 μm (Fig. 9). These samples were obtained by 'extrusion' of the segmented image along the 3 direction. To investigate the influence of interfacial properties (ii), a soft intermittent layer, about 1% of muscle fibre volume fraction (≈ 87%) was introduced between ECM and muscle fibres, and furnished with neo-Hookean type material model by setting̃1 = 0.00155 kPa and̃0 = 1.0 in Eq. (9). The layer is rigidly connected to both the ECM and muscle fibres at their respective interfacing surfaces through shared nodes in the discretised muscle sample.

Results
In what follows, the outcome of the iFEM procedure to determine the ECM material parameters is shown. Thereafter, we present the results of the numerical simulations for 90% pre-compressed ECM samples, for 300 μm cubic muscle samples and for 900 μm extruded muscle samples with and without an intermediate layer between muscle fibres and ECM. Finally, if available, these results are compared to experimental data. Note again that two cases of transverse and oblique loading are considered, respectively, due to the missing symmetry in the microstructure of the sample, providing two response curves (cf. Fig. 5).

Determination of ECM parameters by means of iFEM
The optimisation process was terminated when the objective function had reduced to a value of 0.02. Table 2 lists the obtained parameter set opt , which leads to accurate agreement with the corresponding experimental data for UAE in 0 • , 45 • , and 90 • directions (Fig. 10).

Predictions for compressed ECM
The simulated responses of the laterally compacted ECM samples, generated by shrinking one of the lateral dimensions by 90%, is shown in Fig. 11 in comparison with the experimental data from Kohn et al. (2021) and in Fig. 12 together with the new data on the biaxial response of the compacted ECM. Noteworthy, here the experimental threshold forces of 3 mN for UAE (Kohn et al., 2021) and 10 mN for BAE, with the corresponding respective averaged pre-stress of 0.076 kPa and 0.1 kPa at = 1, were included (observed as a small offset at the stress axis). This is because the simulated compaction operation caused a significantly softer initial response, and thus an extension of the toe region. Even a small threshold force can therefore become relevant and substantially change the appearance of the curve. This soft initial region is also expected to occur in the experiments but R. Kuravi et al.

Material
Parameter Value to be 'hidden' because of the force threshold (see Section 4). When including this threshold, the UAE predictions closely match the experimental response. The BAE predictions, however, though predicting the order of the curves right, underpredict the experimental response.

Mechanical response of cubic muscle tissue samples
The numerically estimated effective nominal stress response of 300 μm cubic muscle samples subject to UAE ( = 1.3) and UAC ( = 0.76) along 0 • , 45 • and 90 • directions to the gross muscle fibre direction ( 3 ) are shown in Figs. 13a and b, respectively. From Fig. 13a, it is observed that loading along 0 • direction provides a weaker response than along the 45 • directions up to a stretch of about 1.17, after which the latter evolves as the stiffer direction (2-3 times stiffer at = 1.3). The weakest response is predicted along the 90 • direction. To facilitate a fair comparison between experiments and simulations, cf. Kohn et al. (2021), the stretches of the computed responses (Fig. 13a) were rescaled such that the experimental initial pre-loads of 30 mN in UAE (equivalent nominal stress of 0.76 kPa) and 10 mN for BAE (equivalent nominal stress of 0.05 kPa) were matched at = 1. (Fig. 13c). The comparison reveals that the ascending/descending order of stiffness, i.e. 90 • > 45 • > 0 • is in line with experimental evidence (Kohn et al., 2021;Takaza et al., 2013a), while there is a mismatch in magnitude in that the predicted 0 • direction is stiffer, whereas the computed 90 • direction is weaker.
Underprediction of the magnitude but a match in the order among the anisotropic stresses is also observed for the BAE response (Fig. 13d) of muscle samples along 0 • and 90 • directions. For instance, at = 1.1 the predicted responses reach only about 30% of the mean experimental observations.
The numerically estimated effective nominal stress responses to modes I, II and III of semi-confined compression (Figs. 5g-l) are illustrated in Fig. 14 for a compressive stretch up to = 0.8, respectively. Figs. 14b-d compare the predicted responses with the corresponding experimental data (Böl et al., 2014). While the responses to modes I and II compression come close to the experimental spread, the mode III R. Kuravi et al. Fig. 10. Comparison of the experimental nominal stress response of the ECM  with the corresponding fitted response along (a) longitudinal (0 • , (b) 45 • , and (c) transverse (90 • ) directions for the optimised material parameter set opt (1040 iterations of the algorithm Fig. 7). Shaded regions represent standard deviation. simulations (Fig. 14d) overpredict the mean experimental stress by up to an order of magnitude. This is also reflected in the order of stiffness observed among the loading directions shown in Fig. 14a, where mode III > II > I, in clear contrast to mode II > I > III reported in the experiments (Böl et al., 2014).

Mechanical response of larger muscle sample
Figs. 15a,b, respectively, compare the UAE and UAC effective nominal stress responses of the extruded 900 μm × 900 μm × 20 μm muscle sample with that of the full 3D 300 μm cubic sample for 0 • and 90 • loading directions. The responses of both the samples are similar both qualitatively and quantitatively. More precisely, in UAE (Fig. 15a) the extruded sample exhibits 20% weaker response in 0 • direction, while it is 15%-30% stronger in the 90 • direction when compared to the 300 μm sample at peak stretch ( = 1.3). Similar deductions are made in UAC (Fig. 15b) wherein, the extruded sample is 44% weaker in 0 • direction and 11%-18% stronger in 90 • direction compared to the 300 μm sample response at peak stretch ( = 0.82).   . Note that two cases of 45 • and 90 • loading were considered, respectively (cf. Fig. 5). Subscript refers to results without considering pre-loads. Fig. 15c compares the effective nominal stress response of the 900 μm extruded sample to the 300 μm sample for all three modes of SCC (Figs. 9e-f), wherein it is observed that the response is similar for modes I and II with marginal differences around ≈ 10%, while in mode III the large model reaches nearly twice as high stress at = 0.82. Fig. 16 compares the behaviour of the extruded sample with and without intermittent layer under UAE, UAC and SCC load cases wherein it is observed that the latter, though weaker, closely approximates the former in both order of magnitude and ascending/descending order of stiffness. In UAE/UAC (Figs. 16a, b), they differ, respectively, by about 4%-8% and 2%-4% along 0 • and 90 • directions at peak stretch. In SCC (Fig. 16c), a difference of 2%, 7.5% and 4% is observed at peak stretch in modes I, II and III, respectively.

Discussion
In this paper, we have combined our recently published approach to generate full 3D FE models of skeletal muscle tissue from histological images and novel experimental data on the mechanical behaviour of the collagenous muscle ECM. This combination put us in a position to calibrate the models of individual components and to perform rigorous computational predictions of the composite tissue response. Although there exist several elaborate studies that employ advanced multi-component microstructure-based models formulated either on a constitutive modelling (Spyrou et al., 2017(Spyrou et al., , 2019Bleiler et al., 2019) or FE simulation basis (Sharafi and Blemker, 2010;Virgilio et al., 2015;Sharafi and Blemker, 2011;Marcucci et al., 2017Marcucci et al., , 2019Spyrou et al., 2017Spyrou et al., , 2019, the first time use of the novel ECM data makes the present approach unique in that it does not contain a single parameter that should be fitted from tissue scale experiments. Accordingly, the simulations are unforgiving, and directly reveal misconceptions and erroneous assumptions when composing the tissue from its main constituents.

ECM model and parameter estimation
The comprehensive new data set (Kohn et al., 2021) suggested a remarkably stiff response when loading the decellularised samples in an angle of 45 • with regard to the prior axis of the fibres. We hypothesised that this behaviour could be caused by significant shear across the ECM layers, enforced through this type of loading. Hence, the constitutive model  that governs the mechanical response of the ECM in each integration point of the FE model, was amended by a shear-dependent term ( ) (Eq. (16)) to incorporate energy associated with shear deformation through an invariant (Section 2.3.2). In fact, the previous version of the SEDF  assumed orthotropic material symmetry with two preferential directions, predominantly inspired by the structure of perimysium (Rowe, 1974;Purslow, 1989), and thus ignoring that endomysium exhibits a relatively random distribution (Trotter and Purslow, 1992;Purslow and Trotter, 1994;Purslow, 2008). This was justified with the relative abundance and comparable tensile stiffness of perimysium (Light et al., 1985;Purslow, 2010). However, endomysium, that coordinates force transmission across individual muscle fibres (Sharafi and Blemker, 2011;Huijing, 1999;Purslow, 2008), is expected to offer high shear resistance to facilitate growth and repair of sarcomeres (Purslow, 2010) and load transfer in intrafascicularly terminating fibres (Sharafi and Blemker, 2010;Purslow, 2002). On the other hand, perimysium is R. Kuravi et al. Fig. 14. Numerically simulated SCC response of 300 μm muscle sample in comparison with experiments. (a)-(c) Comparison of mode I, mode II, and mode III type SCC response with the corresponding experimental data. (d) SCC response of the sample in modes I, II and III. Shaded regions represent standard deviation about their respective experimental mean response. Note that for each mode two cases, respectively, were considered (cf. Fig. 5).  Note that two cases of 90 • loading and SCC were considered, respectively (cf. Fig. 5).
hypothesised to be relatively shear compliant to accommodate fascicle scale shear deformation (Purslow, 1999;Sharafi and Blemker, 2010;Purslow, 2010). The shear term might therefore be interpreted as a means to account for the endomysial properties, and remedy the compliant shear behaviour of the previous model. Clearly, for future developments it would be desirable to model endo-and perimysium with their individual properties. However, as long as there is no specific experimental data on the mechanical behaviour of these components, their separation in the computational model might remain of little use.
Both the experimental and simulated stress-strain curves represent combinations of structural and material responses of the heterogeneous collagenous samples, pervaded by the channels that remain of the muscle cells. These characteristics ruled out the use of an analytical solution to fit the experimental data on the direction-dependent mechanical properties (Kohn et al., 2021). Instead the parameters of the revised constitutive model for ECM (Eq. (16)) were determined through an iFEM procedure (Fig. 7), based on 'virtually decellularised' samples whose direction-dependent mechanical properties were fitted to the experimental data set (Kohn et al., 2021). To avoid the numerical effort associated with a complex self-contact problem, the FE model of the ECM was constructed from the whole domain  instead of the isolated ECM portion, that contains pores that could collapse under load. Instead, the muscle fibres were included but furnished with approximately 100 times weaker material properties. Altogether, the parameter identification proved to work stably, and accurate agreement between the simulations and the experimental data could be achieved (Fig. 10).

Predictions with the calibrated models
The calibrated model was used for predictive simulations. At first we predicted the response of samples that were laterally compacted by 90% before the application of a tensile load (Section 2.4.2). Lateral compaction was mimicked in a simplified procedure, by a purely geometric scaling of the ECM samples along one transverse direction, such that the samples were reduced to 10% of their width. This operation changes the local orientation of the ECM structures, i.e. the vectors , ′ , and in that it generates a nearly planar structure. However, it does not account for other alterations of the structure that arise due to interactions between the collagenous components that have now been brought in close contact. Lateral pre-compression was experimentally observed to have little effect on the order of increasing stiffness (45 • > 90 • ≳ 0 • for UAE and 90 • > 0 • in BAE) which is aptly predicted by the numerical simulations (Figs. 11 and 12). The quantitative comparison between the predicted response and the experiments reveals the tremendous influence of the threshold force used to define the reference state during the tests. When neglecting this apparently small force, the predictions for both UAE (Fig. 11) and BAE (Fig. 12) load cases throughout underestimate the experimental data. However, when rescaling the curves such that the threshold force is reached by the model at a tensile stretch = 1, excellent agreement is observed for the UAE tests, and the predictions for the BAE tests are improved.
The effective nominal stress response of virtual muscle samples to various load cases (Fig. 5) was estimated by subjecting the calibrated FE models to the corresponding BVPs imposed through homogeneous boundary conditions. A comparison of the virtual samples with their experimental counterparts in terms of the response in UAE (Fig. 13c), BAE (Fig. 13d), and SCC (Fig. 14) reveals that the predictions generally match the order of magnitude of the experiments, and that the model predictions on the ascending/descending order of stress are also in line with experiments in UAE (Fig. 13c), BAE (Fig. 13d), and modes I and II of SCC (Figs. 14b,c). Yet, except for UAE in 45 • and BAE along 90 • direction, all simulations generally fail to fall within one standard deviation of the experimental responses. Furthermore a qualitative mismatch is observed in the numerically estimated SCC response in mode III (Fig. 14d) for < 0.9, wherein the predicted stresses are not only higher than experimental observations by an order of magnitude but also stiffer than mode I and mode II stresses.

Missing information and limitations
Although it would probably be straightforward to adapt the model parameters post-hoc in order to achieve better agreement, this was not the goal of the present study. Rather, the detailed virtual experimental set-up was employed to shed light on the potentially missing components and information of multiscale modelling approaches to skeletal muscle tissue. This analysis thus not only points to limitations of the present study but also on general limitations in the available data and experimental procedures.

The strong effect of tare loads and reference states
Small initial loads used to define the reference state are common practice in soft tissue testing. Often they are used as tare loads and even neglected in the representation. While these small forces are needed, e.g. to align the specimen on the one hand, and to reach a reliable range of the force sensors on the other hand, they generate some uncertainty since a substantial part of the 'toe region' (Fung, 1981, Ch. 7) may be lost, due to the strongly non-linear -shaped stress-stretch response exhibited by soft tissues in general (cf. Gao and Desai, 2010;Buerzle and Mazza, 2013). A possibility to overcome this problem is the consideration of the pre-force in the representation of the simulations (e.g. Ehret et al., 2017), i.e. representing the simulated response with respect to the same pre-load, which allows a fair comparison with the corresponding experiments. Therefore, we have chosen this representation in all predictions . Figs. 11 and 13 reveal the strong effect of considering the applied pre-load in that it not only affects the magnitude of the stress values but also brings the 0 • , 45 • and 90 • curves closer together. Noteworthy, the effect of the reference stress was neglected in the iFEM procedure to fit the model to experimental data of uncompressed ECM. This represents a small limitation of the current approach since the effect is also present in these data, albeit less pronounced due to the stiffer initial response of the uncompressed samples.
Another question related to the reference state of a composite material is whether this corresponds to the reference states of its components, cf. Kohn et al. (2021). Theoretically, the components may be residually strained so that their separation leads to different 'stressfree' configurations. But even if not, the pre-load together with the different material characteristics of the components will unavoidably lead to different apparent reference states.

Size effects
Experimental data on muscle tissue behaviour are obtained from experiments performed on samples of typically few millimetres lengths scales (Böl et al., 2014;Van Loocke et al., 2006;Takaza et al., 2013b), that contain tens of fascicles and hence display a relatively homogeneous distribution of ECM, in particular perimysium bands, across the sample cross section in comparison to the 300 μm cubic sample utilised in this work for numerical simulations. To investigate to what extent the quantity and uneven distribution of ECM inherent to the small virtual muscle samples could influence our predictions of muscle tissue behaviour, numerical simulations were performed on a 900 μm × 900 μm × 20 μm extruded sample that encompasses several fascicles. Our previous studies  had shown that despite a generally strong influence of 3D microstructure on the response, the nominal stresses along the loading directions in UAE and UAC were only moderately affected. The extruded sample, for which the thickness becomes irrelevant in FE computations due to symmetry, was hence deemed a reasonable compromise between computationally expensive large full 3D samples and the 300 μm cubic sample.
While clear differences are noticeable between the stresses predicted by the small and large samples in UAE/UAC along 0 • and 90 • directions (Figs. 15a, b) as well as mode I and II SCC, and even a very strong effect in mode III SCC (Fig. 15c), these results indicate that the size of the modelled tissue domain is not the critical factor in explaining the deviations from the experimental responses in Figs. 13 and 14.
Notably, all BVPs in this study were implemented through affine displacement boundary conditions, which can exhibit differences from the state of deformation observed at the boundary in corresponding experiments, in particular for composite tissues. It is known that the effect of the boundary conditions decreases with an increasing computational domain according to St. Venant's principle (Murakami, 2016). We therefore emphasise that the structural effect of considering a larger tissue region may be overlaid by the computational effect of considering larger samples.

Interface effects
Experimental observations and numerical studies (Sharafi and Blemker, 2011) suggest that shear based deformations occur at the fascicle boundaries, i.e. in perimysium and to a lesser degree within the fascicle (Purslow, 1999(Purslow, , 2002, and document a complex linkage between ECM and muscle fibres (Csapo et al., 2020) as well as between intra-ECM layers (Passerieux et al., 2006;Purslow, 2008). For instance, endomysium is reported to form a continuous shared layer between muscle fibres while more discrete 'junction' like connections are observed at the muscle fibre, endomysium and perimysium interface. Such complex manifestations cannot be explored with the continuum model framework used for the ECM. Nevertheless, to obtain qualitative insights into the effect of a potentially compliant interface between muscle fibres and ECM, numerical simulations were performed with the 900 μm extruded sample including a fairly compliant but rigidly bonded hyperelastic intermittent layer. Although the effect may be amplified in full 3D models that favour shear deformations, the marginal effect on the stress (≈ 10%) in all UAE and SCC load cases (Figs. 16ac) delivers a first indication that the interface properties are of minor importance for the predicted stresses in this work.

Non-collagenous ECM, and ECM destruction
Previous studies (Purslow and Trotter, 1994;Passerieux et al., 2007) suggest that the treatment of muscle tissue samples with NaOH as used in Kohn et al. (2021) to isolate ECM not only digests the muscle cell structure but also depletes the ground matrix of the ECM composed of proteoglycans and glycosaminoglycans (PGs/GAGs), or may even cause damage to the fibrous collagen network, depending on concentration and temperature of the solution. The presence of PGs/GAGs could play a particular role in the volumetric response of the material and the low strain regime (see, e.g., Ehret et al., 2017). Moreover, they may play a vital role in collagen cross-linking and maintaining structural organisation (Gillies and Lieber, 2011). In effect, there might be a non-ignorable difference between ECM containing PGs/GAGs and the collagenous ECM alone. Consequently, deviations may exist between the true, i.e. complete ECM behaviour and the one used to parametrise the constitutive model of ECM, particularly with regard to the isotropic contribution to the ECM response.
Partial destruction of the ECM may also explain the early onset of inelastic behaviour of the decellularised muscle samples seen in Kohn et al. (2021). Therein the stress-stretch data to which the ECM model was fitted, is represented within an 'elastic range', within which inelastic behaviour and damage are not expected, and which was smallest for 90 • testing among the uncompressed samples. The stretch that marks the end of this region is given by = 1.2 ± 0.1 (Kohn et al., 2021). To generally avoid exceeding this limit, we restricted the simulations to a smaller stretch range (Figs. 13c,d). However, due to the complex 3D structure the local ECM strains in the muscle FE model may exceed the local strains in the decellularised muscle FE model used for material parameter identification at some locations; whenever this is the case, the ECM model is used beyond the calibrated regime.

Incomplete data
Even if a material model adequately incorporates all relevant mechanisms in a material, the validity of the corresponding material parameters is greatly influenced by the completeness of experimental data used for their identification, e.g. with regard to anisotropy, to the behaviour in tension and compression regimes, and to the tissue response in volumetric and isochoric deformations. Although data from multi-directional loading (0 • , 45 • , 90 • ) was used to obtain the optimal parameter set ( ) of the ECM, the calibration was based on only tensile data and, moreover, no information on the volumetric properties was available. The latter seem of particular importance since the heterogeneity of the sample might favour spatially varying volume changes, so that the volumetric properties of the constituents become crucial even in the case where the muscle sample is little constrained at the tissue scale.

Summary and conclusions
In this work, a recently proposed framework to establish histologybased full 3D models of skeletal muscle tissue was used, adapted and calibrated by recent sets of experimental data on single muscle fibre behaviour and the anisotropic characteristics of decellularised muscle tissue, containing the collagenous ECM. The rigorous predictions of tissue behaviour by this component-level calibrated model capture typical mechanical features of muscle tissue such as anisotropy and tension-compression asymmetry, and also provide reasonable order of magnitude estimates for most load cases. Also, when compared to the mean responses of samples from the same muscle type and species as used for calibration, the predictions fall into the right range, however, accurate agreement is not achieved. While one would need a larger number of virtual samples together with statistics to ensure that the observed deviations are not just a result of the existing variation between the tissue samples used for discretisation, calibration and validation, the deviations point at missing information and ingredients of the bottom-up strategy. Since there are no parameters defined at the tissue scale that could remedy the mismatch, such shortcomings are forthrightly revealed by the present model but they may affect all microstructure-based models. The investigations presented herein indicate that three aspects have a particularly strong effect and should be addressed in future work: The reference states of the single components alone and within the tissue need be clarified, in particular in conjunction with the use of experimental tare loads that characterise the response threshold of the force sensors. The contribution and mechanical characteristics of non-collagenous ECM components, that were most probably lost in the current data set on ECM behaviour, must be quantified and included. Information on the volumetric behaviour of ECM and single fibres is required and, potentially, such data need to be inferred from in-situ testing combined with microscopy, since the separation and preparation of the single components may affect these properties in particular.

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.