2cm2cm2cm2cm

# A three-dimensional model of skeletal muscle tissues††thanks: This work was
partially supported by CONICYT-Chile, the Pacific Institute for the
Mathematical Sciences, and NSERC-Canada through the Discovery program.

Javier A. Almonacid33footnotemark: 3 Sebastián A. Domínguez-Rivera Ryan N.
Konno33footnotemark: 3
Nilima Nigam Stephanie A. Ross Cassidy Tam33footnotemark: 3 James M. Wakeling
Department of Mathematics and Statistics, University of Saskatchewan, Canada
<EMAIL_ADDRESS>Department of Mathematics, Simon Fraser University,
Canada (javiera<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>cassidyt@sfu.ca).Corresponding author<EMAIL_ADDRESS>of
Physical Therapy, University of British Columbia, Canada
(stephanie.ross@ubc.ca)Department of Biomedical Physiology and Kinesiology,
Simon Fraser University, Canada<EMAIL_ADDRESS>
###### Abstract

Skeletal muscles are living tissues that can undergo large deformations in
short periods of time and that can be activated to produce force. In this
paper we use the principles of continuum mechanics to propose a dynamic, fully
non-linear, and three-dimensional model to describe the deformation of these
tissues. We model muscles as a fibre-reinforced composite and transversely
isotropic material. We introduce a flexible computational framework to
approximate the deformations of skeletal muscle to provide new insights into
the underlying mechanics of these tissues. The model parameters and mechanical
properties are obtained through experimental data and can be specified
locally. A semi-implicit in time, conforming finite element in space scheme is
used to approximate the solutions to the governing nonlinear dynamic model. We
provide a series of numerical experiments demonstrating the application of
this framework to relevant problems in biomechanics, and also discuss
questions around model validation.

Keywords: skeletal muscle, non-linear deformation, continuum mechanics,
hyperelasticity, dynamic contraction, isometric contraction, finite element
method, time stepping method

MSC: 74B20, 74L15, 65M60

##  1 Introduction

Why do animals move the way they do? Do the muscles of mice and giraffes and
frogs behave identically? Is the force produced by a muscle (scaled by size)
uniform across species? What transpires during neuromuscular disorders, and of
the many changes which we see – which are the most significant? Such questions
continue to motivate large numbers of researchers, across the fields of
biology, physiology, biomechanics, materials science, and medicine. A full
understanding of the mechanical role of muscle architecture and the properties
of its tissues remains an open and challenging scientific problem. This
article presents some recent work involving the mathematical modeling and
simulation of skeletal muscle mechanics. Of course, given the complexity of
animal locomotion and muscle behavior, numerous important questions are not
addressed in this work. We hope this article sparks the interest of readers
and encourages the involvement of more mathematical research into this
fascinating scientific area.

Skeletal muscles are highly complex tissues which play a crucial structural
and mechanical role in living beings. Unlike non-living materials, muscles can
be deformed both actively or passively. Active deformations are experienced
when internal electro-chemical processes and changes in the shape of proteins
are triggered by neuronal signals (see, e.g. [51, 88]). Passive deformations
describe the non-linear mechanical response of muscle tissue in the presence
of external forces [25]. Both voluntary and involuntary movements are possible
within skeletal muscle; the cells within skeletal muscle are long and multi-
nucleated, and respond to external neuronal stimuli. In contrast, cardiac
muscle (which has also been extensively studied within the mathematics
community) consists of mono-nucleated cells which contract on their own
intrinsic rhythm; see, for instance [46].

Given the complexity of skeletal muscle, teasing apart the role played by
shape (the length, cross-sectional area, angle of attachment to joints),
architecture (the arrangement of constitutive tissues), pennation (the angle
at which muscle fibres are arranged), the arrangement of fast and slow-twitch
fibres, which fibres are enervated and the mechanical properties of the
individual tissues (including their their stiffness, compressibility, density)
is notoriously challenging in experiments. Several of these change as we age –
there is increased fat infiltration in our muscle, we lose muscle mass,
connective tissues become stiffer. Which of these contributes more
significantly to loss of muscle force output, or are these all equally
important effects? Likewise, neuromuscular diseases lead to changes in
architecture, shape, muscle properties and overall muscle mass. It is near-
impossible to vary these properties within a given test subject experimentally
to determine the differential impact of these changes. The goal of therapeutic
treatments is to address the most significant factors impacting muscle
function. At the same time, advances in medical imaging provide us with the
possibility of generating patient-specific computational meshes on which we
can examine the role of these factors. The framework we present in this paper
has been used to study some of these questions, and in this paper we describe
this in detail. One of the studies included here is on integrating MRI-based
data to study muscle contraction; another is on furthering our understanding
of cerebral palsy. We believe that advances in experimental science, medical
imaging and computational power can be powerfully used in tandem to further
our understanding of muscle mechanics.

One-dimensional models of skeletal muscle have been developed to achieve a
deeper understanding of skeletal muscle (e.g. [25, 96]). These models,
generally referred to as Hill-type muscle models, have provided a substantial
insight into the overall behavior of skeletal muscle. This includes insight in
the basics of muscle mechanics [96], the role of muscle in the musculoskeletal
system [14], and the impact of muscular diseases on force output [83]. While
these models are simple and reasonably accurate in some situations
(particularly when considered as part of a full musculoskeletal study), they
fail to account for many features of whole muscle dynamics and architecture.
Mass and inertial effects are not considered in most Hill-type models, and
indeed may be negligible if considering single-fibre experiments. It is
frequently assumed that larger muscles behave as scaled-up single fibres (see,
e.g. [96]). However, in larger muscle the impact of mass cannot be neglected,
especially when the fibres are sub-maximally activated (see eg. [26, 69, 67]).
In addition, the three-dimensional architecture of muscle cannot be
represented by Hill’s model alone. Muscles are nearly incompressible tissues
[7], and thus, when muscle deforms, parts of the surface of the muscle deform
inwards while other parts deform outwards, even during isometric contractions
[60, 71, 88]. In this paper an isometric contraction is defined as an elastic
deformation of the tissues in which two opposing ends of the muscle are fully
clamped. In the physiology literature however, isometric contraction refers to
the combined length of the muscle and tendon being invariant. Additionally,
the deformation also depends on other tissues not explicitly considered in
many one-dimensional models; aponeurosis, base material (comprising
extracellular materials and fluids), and fat all affect the dynamics and
mechanics of muscles, see e.g. [58, 68, 11]. Therefore, to be able to better
capture and understand the behavior of muscle deformation in three dimensions,
a more involved model is needed.

To address the aforementioned concerns and fully capture the mechanics behind
skeletal muscle, three-dimensional models have been developed. Many of these
models are phenomenological continuum models that are developed using whole
muscle properties (e.g. [9, 21, 23, 22, 32, 38, 37, 45, 52, 98, 57, 73, 70,
88, 89, 90, 95]). We note in particular [9] (a quasi-static approach), which
has lead to an explosion of three-dimensional muscle models coupled with
finite element simulations. These models are particularly useful as they are
able to capture the architecture and bulging of muscle. Furthermore, these
models have been developed in both quasi-static and dynamic formulations.
While these models have been able to accurately capture the overall mechanics
of skeletal muscle, they are computationally expensive. This makes them
difficult to implement in larger musculoskeletal models such as OpenSim [14].
Recently, authors in [18] developed a fibre based model, which has a lower
computational cost than the typical continuum models, for the purposes of
implementing more complex simulations with multiple muscles. While these
models focus on the mechanics of the whole muscle, to reach a deeper
understanding of how muscle functions, the influence of the microstructure
needs to be considered.

The length scales involved in muscle mechanics range from the $\mu$m for
sarcomeres to the size of the gastrocnemius in larger mammals (in the order of
meters) it is computationally expensive to model the behavior of large muscles
by considering the individual behaviors of each filament. Micromechanical
models have been developed to capture the influence of features, such as the
connections of the cross-bridges, on the macroscopic scale [10, 36,
62].Another method of capturing microscopic mechanics on a macroscopic scale
is the use of homogenization methods (e.g. [8, 33, 80]). These models are
utilized to capture the influence of the microstructure on the macroscopic
scale in a computationally efficient way. Additionally, some continuum models
have been developed to look at force transmission from microscopic features to
the macroscopic structures [76].

Skeletal muscle models have also been developed to specifically address the
mechanics of impaired muscle tissue, such as that affected by disease, disuse,
or aging. For example, authors in [81] developed a chemo-mechanical model of
muscle tissue that captures the influence of muscular dystrophy, and the
implications of the disease on the fibre level have been investigated in [85].
Continuum models have yielded qualitatively similar results around regional
deformations as experimental data of maximal contractions at constant muscle
velocity, as well as during fixed-end contractions (e.g. [68]). Further, these
models can be used to understand the influence of transverse compression [71]
and fat infiltration [58] on muscle force.

The continuum models of muscle mentioned before focus on quasi-static
deformation, and do not include any velocity-dependent effects; however, the
velocity-dependent creep and relaxation have been found to be important in
muscle mechanics [25]. Viscoelastic models of muscle have been developed to
capture these effects. [84] developed a viscoelastic model to look at the
time-dependent properties of muscle in passive muscle tissue. More recently,
[1] developed a passive viscoelastic model using experimental data for the
cyclic lengthening tests. These models do not, however, look at the active
mechanics of muscle tissue using a viscoelastic model. Other dynamic effects
that need to be considered are the effects of mass [69, 65, 68]. Muscle
deformation is inherently dynamic, and so in this paper we describe a fully
dynamic model of skeletal muscle that will capture the velocity-dependent
contributions to muscle mechanics in three dimensions. Our approach is
focussed on the meso-scale. With this model we are able to capture the
influence of different material types found in the muscle-tendon unit (MTU):
muscle, aponeurosis, fat, and tendon. We have performed computational
investigations involving the first three of these both in physiological and
pathological conditions (see e.g. [66, 58]), while the inclusion of tendon is
part of our future work. Each of these materials possess different mechanical
and physiological properties, and understanding the influence of the
individual materials is critical when considering diseased muscle tissue [34].

Most of the current knowledge on muscle mechanics comes from experiments
performed in laboratories, and generally on mammals. An important part of our
work is to be able to fit the models that describe the mechanics of each of
the tissues in skeletal muscle to experimental data, generally available in
the form of strain-stress data points. These data are collected under
different experimental conditions, and from different species. As a
consequence, we restricted our attention to parameters in the model which were
significant, and have tested the sensitivity of computations to variations in
them. Validation of the overall model is key, particularly to avoid over-
fitting. This is later discussed in more depth in the paper.

In order to properly capture many different biomechanical phenomena, we
designed our numerical framework to be modifiable to a variety of different
experimental settings. An important design principle was to accurately capture
the dynamics, while still being able to model a quasi-static simulation of
muscle contraction. in silica experiments have been performed with the use of
our framework for isometric contractions to study the energetics of muscle
contractions [88], and of compression loads on muscles [71]. In the dynamic
case, the framework has been used to simulate cyclic contraction [64, 68].
With this framework it is possible to specify different mechanical properties
of various tissues (muscle, aponeurosis, fat, base material and tendon),
localized fibre orientation vectors, regionalized activation of muscle fibres,
as well as importing custom-made (MRI-derived) meshes.

The rest of this paper is organized as follows. In section 2 we provide an
overview of muscle physiology, non-linear elasticity, and introduce the system
of equations that we use to describe the large deformation of muscles. We also
discuss the role of muscle mass. In section 3 we describe the methods used to
determine the discretization strategy used for the model equations. A
discussion of parameter fitting to data is presented in section 4. Finally,
some results to validate our framework and demonstrate its biomechanical and
clinical significance are presented in section 5.

##  2 Muscle physiology, mechanics, and modelling

In this section we first overview skeletal muscle physiology and include a
description of the classical Hill-type muscle model [25]. Then we formulate a
three-dimensional dynamic skeletal muscle model, as well as its quasi-static
counterpart.

###  2.1 Skeletal muscle physiology

While skeletal muscle structure varies from micrometers to the scale of
meters, the modelling approach in this study will be focused on a macroscopic
description of muscle and considers muscle fibres as the smallest unit. Muscle
fibres are long, cylindrical cells composed of contractile unites called
sarcomeres arranged in parallel; the sarcomere is the functional unit of
muscle fibres, which contract to produce force. The sarcomeres are on the
scale of micrometers (typically about $2.2\mu$m in length [13]). They consist
of thick myosin filaments interdigitated with thinner actin filaments, which
are connected together via cross-bridges. The force production by the
sarcomere can be characterized using the cross-bridge theory developed by A.F
Huxley [29] and H.E. Huxley[30]. Between each sarcomere is the Z-disk, which
connects to the myosin through a protein called titin. These proteins are
largely responsible for any passive response from the muscle fibres (see e.g.
[24]).

Thousands of muscle fibres are embedded in a matrix of connective tissue
(ECM), which is separated into distinct layers: the epimysium, perimysium, and
endomysium (Figure 1). Each of these layers have a complex structure due to
the collagen fibres; in particular, the collagen fibres have been observed to
have an orientation not necessarily aligned with the muscle fibres [56], which
gives the ECM an anisotropic response. This gives muscle tissue overall a very
complex mechanical response. The fibres render the overall tissue transversely
isotropic, while the ECM introduces its own anisotropy. Typically, a
simplification can be made to model muscle as overall transversly anisotropic
[51, 88, 97] with the anisotropy in the direction of the muscle fibres. The
tendon and aponeurosis play an important role in the performance and
physiology of muscles, as they connect the muscle to the skeletal system.
These tissues are mainly composed of collagen fibres within a matrix of
connective tissue. In tendons, the collagen fibres are arranged in parallel,
running from the bone to the muscle, while aponeuroses are thin sheets of
tendon-like tissue where the collagen fibres are distributed in a more random
fashion.

Figure 1: Structure of a skeletal muscle.

###  2.2 One-dimensional Hill-type muscle model

Most of our current understanding of muscle mechanics comes from measurements
in sarcomeres [48], single fibres, small muscles (such as rat muscle), and in
silica experiments of one-dimensional non-linear models. The latter are based
on the Hill-type muscle model, introduced by A.V. Hill in [25]. This model
provides reasonable estimates of muscle forces during slow contractions [15]
with maximal activation of the muscle fibres [54, 87, 43, 74]. Such models
assume that muscles behave as one-dimensional non-linear damped spring
systems, and that both mass and inertial effects are negligible; in
particular, it assumes that the total force from the muscle is given by the
following equation

$\displaystyle
F_{mus}(\lambda,\epsilon):=F_{0}(a(t)\hat{F}_{act}(\lambda,\epsilon)+\hat{F}_{pas}(\lambda)),$

where $\hat{F}_{act}$ is the muscle force due to the activation of the muscle
fibres (contractile element force), $\hat{F}_{pas}$ is the non-linear elastic
force (passive elastic element), usually referred to as the passive force,
$F_{0}$ is the maximum isometric force of muscles, $\lambda$ is the one-
dimensional stretch, and $\epsilon$ is the rate of change of the stretch over
time, commonly referred to as the strain rate. The function $a=a(t)\in[0,1]$
represents the activation level of the muscle fibres over time. The
definitions of the active and passive forces are normalized to the maximum
isometric force $F_{0}$. We comment that in the physiological community and in
the context of one-dimensional models, the strain rate $\epsilon$ is generally
referred to as the muscle velocity.

A key assumption in Hill’s model is that these forces depend on both the
length and velocity at which the fibres are being stretched. The form of the
active forces is given by

$\displaystyle\hat{F}_{act}(\lambda,\epsilon)=\hat{F}_{len}(\lambda)\hat{F}_{vel}(\epsilon),$
(1)

where $\hat{F}_{len}$ are the forces generated from length changes in the
muscle and $\hat{F}_{vel}$ are the forces due to the changes in the strain
rate. The force $\hat{F}_{vel}$ satisfies $\hat{F}_{vel}(0)=1$. The forces
$\hat{F}_{len}$ and $\hat{F}_{vel}$ are defined as non-linear functions of the
stretch $\lambda$ and the strain rate $\epsilon$ respectively, as seen in
Figure 2. These are commonly called force-length and force-velocity
relationships. Meanwhile, the passive forces $\hat{F}_{pas}$ act to resist
fibres from stretching to longer lengths. These play an important role in
muscle contraction as non-linear springs: when muscle fibres are lengthened,
non-linear passive forces are produced to shorten the fibre back. However,
when fibres contract to shorter lengths, experimental data suggest that the
passive forces in this regime are negligible, at least for small deformations
(Figure 2); see, for instance [68].

###  2.3 A three-dimensional hyperelastic model for skeletal muscle

Even though one-dimensional models are able to explain the behavior of whole-
muscle mechanics in some specific cases, there are many other cases in which
full three-dimensional mechanics is needed to better explain how muscles
behave. This is the main focus of the present section.

####  2.3.1 Basics of continuum mechanics

Let us start by recalling some basic concepts in the theory of non-linear
elasticity. This overview is mainly based on [27, 50, 86]. Suppose that an
elastic body initially occupies a bounded region $V_{0}\subset\mathbb{R}^{3}$.
We consider, for fixed $t>0$, that there exists a diffeomorphism
$\phi(\cdot,t):V_{0}\to V$ which maps every point in $V_{0}$ onto a bounded
region $V\subset\mathbb{R}^{3}$, that is the range of $\phi(\cdot,t)$ is
exactly $V$. The region $V_{0}$ is referred to as the initial (or reference)
configuration on the elastic body while $V$ represents the current
configuration of such a body. We assume that this map is such that
$\phi(\cdot,0)$ is the identity function. Thus, a point ${\bf x}\in V$ can be
written as ${\bf x}=\phi({\bf X},t)$, for some ${\bf X}\in V_{0}$. We define,
at time $t>0$, the material displacement of a point ${\bf X}\in V_{0}$ and the
spatial displacement of a point ${\bf x}\in V$ as

$\displaystyle{\bf U}({\bf X},t):={\bf x}-{\bf X}=\phi({\bf X},t)-{\bf
X},\quad{\bf u}({\bf x},t):={\bf x}-{\bf X}=x-\phi^{-1}({\bf x},t).$

For points ${\bf X}\in V_{0}$ and ${\bf x}\in V$, the material velocity and
the spatial velocity are given respectively by

$\displaystyle{\bf V}({\bf X},t):=\frac{\partial\phi({\bf X},t)}{\partial
t}=\dfrac{\mathrm{D}{\bf U}}{\mathrm{D}t}\equiv\dot{{\bf U}},\quad{\bf v}({\bf
x},t):={\bf V}(\phi^{-1}({\bf x},t),t)=\dfrac{\partial{\bf u}}{\partial t}.$

where $\frac{\mathrm{D}}{\mathrm{D}t}\equiv\dot{(\cdot)}$ denotes the material
time derivative of a material field and $\frac{\partial}{\partial t}$ denotes
the spatial time derivative of a spatial field. We also define the deformation
gradient tensor and its determinant (also known as the Jacobian determinant)
at a given material point to be

$\displaystyle{\bf F}({\bf X},t):={\bf I}+\nabla_{0}{\bf U}({\bf X},t),\quad
J({\bf X},t):=\det{\bf F}({\bf X},t),\quad{\bf X}\in V_{0},\,\,t>0.$

As for strain measures, we consider the right and left Cauchy-Green tensors,
respectively given by

${\bf C}={\bf F}^{\texttt{t}}{\bf F},\quad{\bf b}={\bf F}{\bf
F}^{\texttt{t}}.$

We also introduce the material velocity gradient and the spatial velocity
gradient:

$\dot{{\bf F}}=\dfrac{\partial}{\partial t}\left(\dfrac{\partial\phi({\bf
X},t)}{\partial{\bf X}}\right)=\dfrac{\partial{\bf V}({\bf X},t)}{\partial{\bf
X}},\quad\l:=\dfrac{\partial{\bf v}({\bf x},t)}{\partial{\bf X}}=\dot{{\bf
F}}{\bf F}^{-1}.$

The symmetric part of $\l$ corresponds to what is known as the rate of strain
tensor ${\bf d}$, for which the following holds:

${\bf d}:=\dfrac{1}{2}\left(\l+\l^{\texttt{t}}\right),\quad\dot{{\bf C}}=2{\bf
F}^{\texttt{t}}{\bf d}{\bf F}.$

Finally, we define the stretch $\lambda$ and strain rate $\epsilon$ of fibres
with respect to the reference configuration as

$\lambda=\|{\bf F}{\bf a}_{0}\|,\quad\epsilon:=\dot{\lambda}=\dfrac{({\bf
F}{\bf a}_{0})^{\texttt{t}}\,{\bf d}\,{\bf F}{\bf a}_{0}}{\lambda},$ (2)

where ${\bf a}_{0}\in\mathbb{R}^{3}$ is a unit vector representing the initial
orientation of a fibre at a point ${\bf X}\in V_{0}$ and $\|\cdot\|$ is the
usual 2-norm of vectors in $\mathbb{R}^{3}$.

####  2.3.2 Multiplicative decomposition of the deformation tensor

In the case of compressible materials (even by small amounts, i.e. nearly
incompressible), it is customary to split the deformation tensor into volume-
changing (volumetric) and volume-preserving (isochoric) parts, i.e.

${\bf F}=J^{1/3}\bar{{\bf F}},\quad\bar{{\bf F}}:=J^{-1/3}{\bf F}.$ (3)

We call $\bar{{\bf F}}$ the modified deformation gradient and we have that
$\det\,\bar{\bf F}=1$ (this is the isochoric part of the deformation). This
decomposition yields modified versions of the right and left Cauchy-Green
tensors:

$\bar{{\bf C}}:=\bar{{\bf F}}^{\texttt{t}}\bar{{\bf F}}=J^{-2/3}{\bf
C},\quad\bar{{\bf b}}:=\bar{{\bf F}}\bar{{\bf F}}^{\texttt{t}}=J^{-2/3}{\bf
b},$ (4)

as well as modified stretch and strain rates:

$\bar{\lambda}:=J^{-1/3}\lambda,\quad\bar{\epsilon}:=\dot{\bar{\lambda}}=\dfrac{(\bar{{\bf
F}}{\bf a}_{0})^{\texttt{t}}(\mathbbm{p}:{\bf d})(\bar{{\bf F}}{\bf
a}_{0})}{\bar{\lambda}},$ (5)

where $\mathbbm{p}:=\mathbb{I}-\frac{1}{3}{\bf I}\otimes{\bf I}$ is the
spatial projection tensor and $\mathbb{I}$ is the fourth-order symmetric
identity tensor, i.e. in index notation
$\mathbb{I}_{abcd}=\frac{1}{2}(\delta_{ac}\delta_{bd}+\delta_{ad}\delta_{bc})$.

Several invariants (and pseudo-invariants) have been defined for the right
Cauchy-Green tensor ${\bf C}$ [91]. In particular, the models to be presented
will make use of the first and fourth modified invariants of $\bar{{\bf C}}$,
with derivatives indicated as follows:

$\bar{I}_{1}={\rm tr}\,\bar{\bf C},\quad\bar{I}_{4}={\bf
a}_{0}^{\texttt{t}}\bar{\bf C}{\bf
a}_{0}=\bar{\lambda}^{2},\quad\dfrac{\partial\bar{I}_{1}}{\partial\bar{\bf
C}}={\bf I},\quad\dfrac{\partial\bar{I}_{4}}{\partial\bar{\bf C}}={\bf
a}_{0}\otimes{\bf a}_{0}.$ (6)

The fourth invariant $\bar{I}_{4}$ will be specially important as it contains
information on the direction of the fibres. Moreover, its rate of change,
$\dot{\bar{I}}_{4}=2\bar{\lambda}\dot{\bar{\lambda}}=2\bar{\lambda}\bar{\epsilon}$,
depends directly on the strain rate of the fibres.

####  2.3.3 A strain-energy function for the muscle-tendon unit

As described previously, our model of skeletal muscle consists of
contributions from the muscle, tendon, fat, base material, and aponeurosis. We
consider muscle, aponeurosis, and tendon as anisotropic fibre-reinforced
composite and nearly incompressible materials [49]. Consequently, the
Helmholtz strain-energy function $\Psi$ of the material (as a whole) can be
described through a splitting into volumetric and isochoric parts:

$\Psi=\Psi_{vol}(J)+\Psi_{iso}(\bar{{\bf C}},{\bf a}_{0}\otimes{\bf
a}_{0};\dot{\bar{{\bf C}}}).$ (7)

The isochoric part of the strain-energy, $\Psi_{iso}$, encloses the properties
of incompressible deformations through the modified right Cauchy-Green tensor
(recall that $\det\bar{\bf C}=1$). Moreover, $\Psi_{iso}$ contains the
contributions of all parts of the MTU:

$\Psi_{iso}=\Psi_{iso,mus}(\bar{I}_{1},\bar{I}_{4};\dot{\bar{I}}_{4})+\Psi_{iso,apo}(\bar{I}_{1},\bar{I}_{4})+\Psi_{iso,ten}(\bar{I}_{1},\bar{I}_{4}).$
(8)

Here, the dependence on $\bar{I}_{4}$ reflects the fact that the tissues
(muscle, aponeurosis, and tendon) are composed mostly of curvilinear fibres
running from one point to another within $V_{0}$. In the current model it is
assumed that these fibres are one-dimensional (in the along-fibre direction).
However, the space between fibres is not empty and is made of a base material
comprising extracellular material, fluid and other constitutive tissues. The
mechanical response of this material include a response from the extracellular
material, and the three-dimensional cellular response. This adds a dependence
on $\bar{I}_{1}$ in the isochoric strain-energy above. Finally, to include
viscoelastic effects in this model, some dependence on $\dot{{\bf C}}$ must be
added. We propose to include this rate of deformation as a parameter in the
isochoric component of the strain-energy function. This can be justified from
a numerical analysis perspective: in the resulting system of PDEs, an explicit
treatment of the velocity in the time discretization means that the tensor
$\dot{\bar{\bf C}}$ (and consequently $\dot{\bar{I}}_{4}$) is completely known
at the current time step. Thus, the stress and elasticity tensors depend on
the current displacement (and possibly other current variables) but not on the
current velocity. In this way, the fully-dynamic problem can be seen as a
sequence of quasi-static deformations, as a pre-computed velocity is used to
drive the system from current time step to the next one.

Each one of the contributions in (8) can be further subdivided to consider the
different materials found in the MTU, although we save precise definitions for
each strain-energy function for section 4. First, for muscle material, the
micro-mechanical effects can be homogenized to include the ECM, cellular
material, and fat response. In this work, we utilize the breakdown of the base
material derived in [33] with the addition of a volume fraction encompassing
the effects of the fatty tissue:

$\begin{split}\Psi_{iso,mus}&=(1-\beta)\Big{\\{}\alpha\Psi_{ecm}(\bar{I}_{1})+(1-\alpha)\Psi_{cell}(\bar{I}_{1})\\\
&\hskip
100.00015pt+\Psi_{fibre,mus}(\bar{I}_{4};\dot{\bar{I}}_{4})\Big{\\}}+\beta\Psi_{fat}(\bar{I}_{1}),\end{split}$
(9)

where $\alpha$ is the volume fraction of the ECM and $\beta$ is the volume
fraction of fat. Then, for the aponeurosis and tendon, we breakdown the
ischoric component into base material and fibre components:

$\Psi_{iso,tis}=\Psi_{fibre,tis}(\bar{I}_{4})+\Psi_{base,tis}(\bar{I}_{1}),\quad
tis\in\\{apo,ten\\}.$ (10)

We note here that, while the general formulation of $\Psi_{fibre,tis}$ has a
dependency on the strain rate $\bar{\epsilon}$ of the fibres (i.e. on
$\dot{\bar{I}}_{4}$ through (5)), in practice only the muscle fibres will have
this dependency.

####  2.3.4 Potential formulation and dynamic equilibrium

The quasi-static part of the problem can be represented by the stationarity of
the following potential function [77]:

$\Pi({\bf U},p,D)=\int_{V_{0}}\Psi_{vol}(D)+p\left(J({\bf
U})-D\right)+\Psi_{iso}(\bar{\bf C})-\int_{V_{0}}{\bf f}_{0}\cdot{\bf U}.$
(11)

A change in volume during the deformation of the tissues implies a change in
internal pressure, which we denote by $p$. The variable $D$ is defined as the
true dilation $J({\bf U})$ in the reference configuration. Thus, the pressure
$p$ plays the role of a Lagrange multiplier.

Notice that (11) is defined using quantities in the reference configuration.
Thus, a dynamic extension of the Euler-Lagrange equations that arise from the
stationarity of this potential, represent a total-Lagrangian description of
the system:

$\displaystyle\dot{\bf U}={\bf V},\quad\rho_{0}\dot{\bf V}={\bf Div}\,{\bf
P}({\bf U},{\bf V},p)+{\bf f}_{0}\quad\text{in }V_{0},$ (12a) $\displaystyle
J({\bf U})-D=0,\quad p-\Psi_{vol}^{\prime}(D)=0\quad\text{in }V_{0}.$ (12b)

where ${\bf Div}(\cdot)$ is the divergence operator in material coordinates
(i.e. $({\bf Div}\,{\bf P})_{A}=\frac{\partial P_{AB}}{\partial X_{B}}$),
$\rho_{0}=\rho_{0}({\bf X})$ is the density of the tissues in the reference
configuration, ${\bf f}_{0}$ is a body force, and ${\bf P}$ is the first
Piola-Kirchhoff (PK1) tensor. The second Piola-Kirchhoff (PK2) tensor, ${\bf
S}$, is related to the PK1 tensor by the identity ${\bf P}={\bf F}{\bf S}$ and
contains the constitutive laws of the material:

${\bf S}({\bf U},{\bf V},p)={\bf S}_{vol}(p)+{\bf S}_{iso}({\bf U};{\bf
V})=pJ{\bf C}^{-1}+J^{-2/3}\mathbb{P}:\bar{{\bf S}}.$ (13)

Here, $\mathbb{P}:=\mathbb{I}-\frac{1}{3}{\bf C}^{-1}\otimes{\bf C}$ is the
fourth-order material projection tensor and $\bar{\bf S}$ is the ficticious
Piola-Kirchhoff stress given by

$\bar{\bf S}=2\dfrac{\partial\Psi_{iso}(\bar{\bf C};\dot{\bar{\bf
C}})}{\partial\bar{\bf
C}}=2\left(\dfrac{\partial\Psi_{iso}(\bar{I}_{1})}{\partial\bar{I}_{1}}\,\dfrac{\partial\bar{I}_{1}}{\partial\bar{\bf
C}}+\dfrac{\partial\Psi_{iso}(\bar{I}_{4};\dot{\bar{I}}_{4})}{\partial\bar{I}_{4}}\,\dfrac{\partial\bar{I}_{4}}{\partial\bar{\bf
C}}\right).$ (14)

To reveal the inertial component of the system, we pull (12) to the current
configuration. This constitutes the Eulerian description of the system:

$\displaystyle\dfrac{\partial{\bf u}}{\partial t}={\bf
v},\quad\dfrac{\partial(\rho{\bf v})}{\partial t}+{\bf div}(\rho{\bf
v}\otimes{\bf v})={\bf div}\,{\bm{\sigma}}({\bf u},{\bf v},p)+{\bf
f}\quad\text{in }V,$ (15a) $\displaystyle J({\bf u})-D=0,\quad
p-\Psi_{vol}^{\prime}(D)=0\quad\text{in }V,$ (15b)

where the body force ${\bf f}=J^{-1}{\bf f}_{0}$. In particular, inertial
forces in the system are represented by the convective term ${\bf
div}(\rho{\bf v}\otimes{\bf v})$. The Cauchy stress tensor ${\bm{\sigma}}$ can
be written in terms of the Kirchhoff stress ${\bm{\tau}}$ as:

$J{\bm{\sigma}}=:{\bm{\tau}}={\bm{\tau}}_{vol}(p)+{\bm{\tau}}_{iso}({\bf
u};{\bf v})=pJ{\bf I}+\mathbbm{p}:\bar{\bm{\tau}},$ (16)

where $\bar{\bm{\tau}}:=\bar{\bf F}\bar{\bf S}\bar{\bf F}^{\texttt{t}}$ is the
ficticious Kirchhoff stress and $\mathbbm{p}$ is the spatial projection tensor
defined as for (5).

We close the system given by (15) (and equivalently (12)) with initial
conditions and mixed boundary conditions. For this, we assume that the
boundary in the current configuration $\partial V$ is divided as $\partial
V=S_{D}\cup S_{N}$, $S_{D}\cap S_{N}=\emptyset$. The first part, $S_{D}$,
corresponds to the part of the boundary for which displacements are set (e.g.
a clamped face or faces that moves with a prescribed displacement). In turn,
$S_{N}$ represents the set of traction-free surfaces. Mathematically we can
write this as

$\begin{split}{\bf u}={\bf u}_{0},\quad{\bf v}={\bf v}_{0},\quad p=p_{0},\quad
D=D_{0}\quad\text{in }V_{0},\\\ {\bf u}={\bf u}_{D}\quad\text{on
}S_{D},\quad{\bm{\sigma}}({\bf u},{\bf v},p){\bf n}=\mathbf{0}\quad\text{on
}S_{N}.\end{split}$ (17)

####  2.3.5 Quasi-static model

In the quasi-static formulation of the model, we assume that the tissue
experiences slow excitations or deformations, on time-scales much slower than
that for instantaneous mechanical equilibration. In this case, we do not need
the dynamic extension made in (12). Therefore, in the Eulerian description,
the system (15) reduces to:

$\displaystyle{\bf div}\,{\bm{\sigma}}({\bf u},p)+{\bf f}({\bf
x},t)=\mathbf{0}\quad\text{in }V,$ (18a) $\displaystyle D-J({\bf u})=0,\quad
p-\Psi_{vol}^{\prime}(D)=0\quad\text{in }V,$ (18b) $\displaystyle{\bf u}={\bf
u}_{D}\quad\text{on }S_{D},\quad{\bm{\sigma}}({\bf u},p){\bf
n}=\mathbf{0}\quad\text{on }S_{N},$ (18c) $\displaystyle{\bf u}={\bf
u}_{0},\quad p=p_{0},\quad D=D_{0}\quad\text{in }V_{0}.$ (18d)

These are precisely the Euler-Lagrange equations for the stationarity of the
potential (11).

In section 3 we introduce the numerical scheme to approximate the solutions of
the non-linear system in Equation 17. A semi-implicit time stepping is applied
at the continuous level to discretize the time variable. The non-linear
equations are solved via the application of a Newton’s iteration. Finally we
discretize the space variables of the resulting linear system with a mixed
FEM.

####  2.3.6 Muscle mass

We highlight a key question motivating the present work: does muscle mass
matter when describing muscle force output and energy? The immediate response
may be affirmative. However, these effects may be negligible in certain
settings; certainly, in the seminal paper concerning skeletal muscle mechanics
[25], inertial effects due to muscle mass are neglected. These effects are
also negligible, to leading order, when dealing with single-fibre experiments,
or slow deformations. It is common in the literature to use a quasi-static
approximation in these cases. However, it is also known that Hill-type models
(neglecting inertial effect) result in erroneous predictions of actual muscle
force in muscles. Recent experimental evidence [26, 67] demonstrates the non-
negligible impact of muscle mass. In the first of these papers, the authors
consider a rat plantaris muscle in which at first only fast-twitch fibres are
activated. Next, only slow-twitch fibres are activated. One may expect that
when both types of fibres are activated that the maximum velocity of
contraction would be slower than for the fast-twitch types; in fact, it was
seen that activating all the fibres lead to an increase in this velocity. The
authors suggest this is because when only one type of fibre is activated, the
internal mass of the others leads to a reduction in net work output. In [67],
in-situ experiments on rat plantaris muscle were used to study the effects of
increasing the mass of muscle on mechanical work. It was seen that the muscles
with added mass were able to perform less mechanical work, suggesting again
the non-negligible role played by muscle mass in mechanics.

A goal in our work has been to understand if - and under what circumstances -
such mass effects need to be included. As discussed earlier, to leading order
mass effects are perhaps negligible if we model the entire musculoskeletal
system (the rigidity of the skeleton plays a key role), or in single muscle
fibres, or at maximal activation. Recent experiments [35] quantify both shear
and sound wave speeds in tissue; elastic shear waves are seen to propagate at
around 5$ms^{-1}$ (similar to the speed of some neuronal activation pulses),
while sound wave speeds are closer to $1.5\times 10^{3}ms^{-1}$. This clearly
shows that elastic waves – and tissue density – cannot be neglected for fast
twitches.

A simple example will illustrate the point. A common experiment in muscle
physiology consists of clamping a muscle fibre of initial length $L$ at one
end, and pulling it passively with a prescribed displacement ${\bf
u}_{D}=g(t){\rm e}_{1}$ at the other. In the quasi-static approximation, all
points along the muscle are assumed to respond instantaneously to this
displacement, and no compressive elastic wave is expected. In the dynamic
model, however, prescribing ${\bf u}_{D}=g(t){\rm e}_{1}$ becomes equivalent
to adding a force proportional to $\frac{X}{L}\ddot{g}(t)$, where $X$ is the
first component of a point in $V_{0}$. If we, for instance, start to pull the
muscle at a constant velocity and then stop doing so, the system experiences
an impulsive force. The response will include an elastic wave propagating
along the length of the muscle. In section 5, we contrast the computed
mechanical response of the quasi-static and dynamic models.

##  3 Discretization

In this section we describe the numerical methods used to simulate the three-
dimensional deformation of skeletal muscles. Since the size of the nonlinear
system we solve is substantial, we have at this juncture avoided schemes which
add to the list of unknowns (by augmenting the system with the unknown
velocity) that are solved at each time-step. We note other time-stepping
schemes are possible, and are the subject of future work.

###  3.1 Eulerian versus Lagrangian descriptions

At a continuous level, the Lagrangian system (12) and the Eulerian system (15)
are equivalent, with the interchange between integrals and time derivatives
being the link between the two descriptions [27]. Therefore, when the time
variable is discretized first, these two descriptions are no longer
equivalent, but because they are descriptions of the same physical process,
the end result is very similar.

To avoid the complications brought by the convective terms in the Eulerian
description, we choose to discretize the Lagrangian description (12). Notice
also that, unlike in the Eulerian description where the domain $V$ changes at
every time step, the Lagrangian discretization only requires a mesh of the
reference configuration $V_{0}$ that is kept fixed throughout the entire
simulation.

###  3.2 Total Lagrangian formulation

To approximate a solution of the non-linear dynamic system in (15), we use a
semi-implicit time stepping for the time derivatives, and a conforming finite
element scheme for the spatial variables. The displacement, dilation, and
pressure are updated with an implicit approach while the velocity is computed
with an explicit update. The velocity is treated explicitly and it is only
updated once the new displacement is computed. Using this time stepping in
(12), we propose the problem: at time step $t_{n}$, find a material
displacement ${\bf u}^{n}$ (we use lowercase letters from now on for a more
pleasant reading), a material pressure $p^{n}$, and a material dilation
$D^{n}$ such that

$\left\\{\begin{gathered}\rho_{0}{\bf u}^{n}-\delta t^{2}{\bf Div}\,{\bf
P}({\bf u}^{n},{\bf v}^{n-1},p^{n})=\rho_{0}{\bf u}^{n-1}+\delta t\rho_{0}{\bf
v}^{n-1}+\delta t^{2}{\bf f}_{0}^{n},\\\ J({\bf u}^{n})-D^{n}=0,\quad
p^{n}-\Psi_{vol}^{\prime}(D^{n})=0,\end{gathered}\right.$

where $\delta t>0$ is the time step size. The first equation comes from
combining the two equations in (12a) using the time discretization ${\bf
v}^{n}=\left({\bf u}^{n}-{\bf u}^{n-1}\right)/\delta t$, which is also used to
compute the material velocity as a post-process at each time step.

Multiplying these equations by a test function $\delta{\bm{\xi}}:=(\delta{\bf
u},\delta p,\delta D)\in\mathcal{X}$, with $\mathcal{X}:={\bf
H}^{1}(V_{0})\times L^{2}(V_{0})\times L^{2}(V_{0})$, and integrating by parts
one obtains the following non-linear variational formulation: given the
previous displacement ${\bf u}^{n-1}$, previous velocity ${\bf v}^{n-1}$,
force ${\bf f}_{0}^{n}$, and a prescribed displacement ${\bf u}_{D}^{n}$, we
are to find ${\bm{\xi}}^{n}:=({\bf u}^{n},p^{n},D^{n})\in\mathcal{X}$ such
that ${\bf u}^{n}={\bf u}^{n}_{D}$ on $S_{D}$ and

$R({\bm{\xi}}^{n},\delta{\bm{\xi}}):=R_{ine}({\bf u}^{n},\delta{\bf
u})+R_{int}({\bm{\xi}}^{n},\delta{\bm{\xi}})-R_{ext}({\bf u},\delta{\bf
u})=0,\quad\forall\,\delta{\bm{\xi}}\in\mathcal{X},$ (19)

that is, for all virtual displacements $\delta{\bf u}$, virtual pressures
$\delta p$, and virtual dilations $\delta D$. Such a formulation can be
understood as a semi-discrete version of the principle of virtual work in the
total Lagrangian formulation. The semi-discrete work due to inertial,
internal, and external forces are given respectively by:

$\displaystyle\begin{split}&R_{ine}({\bf u}^{n},\delta{\bf u})=\delta
t^{-2}\left(\rho_{0}{\bf u}^{n},\delta{\bf u}\right)_{V_{0}}\\\ &\hskip
80.00012pt-\delta t^{-2}\left(\rho_{0}{\bf u}^{n-1},\delta{\bf
u}\right)_{V_{0}}-\delta t^{-1}\left(\rho_{0}{\bf v}^{n-1},\delta{\bf
u}\right)_{V_{0}},\end{split}$ (20a)
$\displaystyle\begin{split}&R_{int}({\bm{\xi}}^{n},\delta{\bm{\xi}})=\left({\bm{\tau}}({\bf
u}^{n},{\bf v}^{n-1},p^{n}),\nabla_{0}(\delta{\bf u}){\bf F}({\bf
u}^{n})^{-1}\right)_{V_{0}}\\\ &\hskip 80.00012pt+\left(J({\bf
u}^{n})-D^{n},\delta
p\right)_{V_{0}}+\left(\Psi_{vol}^{\prime}(D^{n})-p^{n},\delta
D\right)_{V_{0}},\end{split}$ (20b) $\displaystyle\begin{split}&R_{ext}({\bf
u}^{n},\delta{\bf u})=\left({\bf f}_{0}^{n},\delta{\bf
u}\right)_{V_{0}}.\end{split}$ (20c)

The non-linear nature of this problem is handled using Newton’s method. At
each time step, given ${\bf u}^{n-1},{\bf v}^{n-1},{\bf
u}_{k}^{n},p_{k}^{n},D_{k}^{n}$, we want to find an update ${\bf
d}{\bm{\xi}}_{k}:=({\bf d}{\bf u}_{k},dp_{k},dD_{k})$ such that ${\bf d}{\bf
u}_{0}={\bf u}^{n}-{\bf u}^{n-1}$ on $S_{D}$, ${\bf d}{\bf u}_{k}={\bf 0}$ on
$S_{D}$ for $k\geq 1$, and

$\left[\mathcal{D}R({\bm{\xi}}^{n}_{k})({\bf
d}{\bm{\xi}}_{k}),\delta{\bm{\xi}}\right]=-\left[R({\bm{\xi}}_{k}^{n}),\delta{\bm{\xi}}\right],$
(21)

for all virtual increments $\delta{\bm{\xi}}\in\mathcal{X}$, to then update
the solution as ${\bm{\xi}}_{k+1}^{n}={\bm{\xi}}_{k}^{n}+{\bf
d}{\bm{\xi}}_{k}$, for all $k\geq 0$ until convergence. At the current time
step, we set the initial solution for the Newton iteration as the converged
solution from the previous step, i.e. ${\bm{\xi}}_{0}^{n}={\bm{\xi}}^{n-1}$.
The tangent operator, $\mathcal{D}R$, corresponds to the Gateaux derivative of
the operator $R$ defined in (19), and $[\cdot,\cdot]$ denotes the duality
pairing between $\mathcal{X}$ and $\mathcal{X}^{*}$. For more information on
the linearization of this operator and the resulting fourth-order elasticity
tensors, we refer the reader to [27, 93].

###  3.3 Finite element discretization

To discretize the linear operator $\mathcal{D}R({\bm{\xi}}_{k}^{n})$, we use a
standard conforming finite element. Let $\mathcal{T}_{h}$ be a regular mesh of
$\overline{V_{0}}$ made of (isoparametric) brick elements $T$ of diameter
$h_{T}$. The meshsize is simply defined as
$h:=\max\left\\{h_{T}:T\in\mathcal{T}_{h}\right\\}$. We use the same finite
element spaces as those used in [53], that is, ${\bf Q}_{k+1}\times
P_{k}\times P_{k}$. For a non-negative integer $k$ and an element $T$ of
$\mathcal{T}_{h}$, let us consider the $P_{k}(T)$ as the space of all
polynomials of degree at most $k$ defined in $T$, with ${\bf P}_{k}(T)$
denoting the vector version of $P_{k}(T)$. Then

$\displaystyle{\bf Q}_{k+1}$ $\displaystyle:=\left\\{{\bf
w}\in\mathcal{C}(\overline{V_{0}})\,:\,{\bf w}\big{|}_{T}\in{\bf P}_{k+1}(T),\
\forall\,T\in\mathcal{T}_{h}\right\\},$ $\displaystyle P_{k}$
$\displaystyle:=\left\\{q\in L^{2}(V_{0})\,:\,q\big{|}_{T}\in P_{k}(T),\
\forall\,T\in\mathcal{T}_{h}\right\\}.$

Normally, we will use $k=1$, that is, displacements are approximated using
quadratic 27-node Lagrange elements ${\bf Q}_{2}$, while pressures and
dilations are approximated using linear discontinuous elements $P_{1}$ (based
on monomials, 4 degrees of freedom per element). We also note that the element
${\bf Q}_{1}\times P_{0}\times P_{0}$ does not exhibit locking though the
theory is lacking.

##  4 Determination of model parameters and strain-energy description

Here we describe the methods used to determine the parameters in the material
response of the muscle tissues. The full list of parameters and functions are
included in the Supplementary Material.

###  4.1 Volumetric response

In this study we adopt the form of the volumetric stress used in [53, 49] (see
also [27] for more details on volume changing free energy functions), which is
typically used for nearly incompressible biological materials [27]:

$\Psi_{vol}(J)=\dfrac{\kappa}{4}\left(J^{2}-2\ln\,J-1\right),\quad\kappa({\bf
X})=\sum_{tis\in\\{mus,apo,ten\\}}\kappa_{tis}{\bm{\chi}}_{tis}({\bf X}),$

in which $\kappa_{tis}>0$ is the bulk modulus of the corresponding tissue and
${\bm{\chi}}_{tis}$ is the indicator function of the tissue subdomain. In
particular, the muscle bulk modulus is computed as
$\kappa_{mus}=(1-\beta)\left(\alpha\kappa_{ecm}+(1-\alpha)\kappa_{cell}\right)+\beta\kappa_{fat}$.
The values for bulk moduli are determined through ensuring that the near
incompressibility condition is met, ie. the values for the moduli are chosen
so that there is little change in the volume of the muscle. Based on [9, 57],
we choose $\kappa_{ecm}=1\times 10^{6}$Pa. Due to the incompressibility of
water and the fact cells are mainly composed of water we choose the volume
fraction for the cellular material and fat to be
$\kappa_{cell}=\kappa_{fat}=1\times 10^{7}$Pa [33]. It has been shown in [17,
33] that the resulting strain mechanics are not very sensitive to the bulk
moduli, so we mainly choose these values to ensure that the volume change is
small. Due to the increased stiffness of the aponeurosis and tendon materials,
the bulk modulus is chosen to be larger, $\kappa_{apo}=\kappa_{ten}=1\times
10^{8}$Pa.

###  4.2 Isochoric response due to base material properties

The mechanical responses of the isochoric components of the base material for
each of the different tissues (including ECM and cellular components for
muscle) are modeled using Yeoh-type strain-energy functions [94]. In turn, the
intramuscular fat was modeled as a material of Neo-Hookean type. This yields
two types of strain-energies that will be used in this work:

$\Psi_{Yeoh}(\bar{I}_{1})=\sum_{k=0}^{3}c_{k}\left(\bar{I}_{1}-3\right)^{k},\quad\Psi_{Neo}(\bar{I}_{1})=c_{1}\left(\bar{I}_{1}-3\right).$
(22)

where the constants $c_{k}\in\mathbb{R}$, $k=0,\dots,3$, are determined based
on experimental data for each one of the tissues. For the Yeoh model, the
basis of the fitting procedure will consider the incompressible version of the
Cauchy stress response:

${\bm{\sigma}}_{Yeoh}=-p_{Yeoh}{\bf
I}+2\dfrac{\partial\Psi_{Yeoh}(I_{1})}{\partial I_{1}}{\bf b},\quad
p_{Yeoh}:=\dfrac{2}{\lambda}\dfrac{\partial\Psi_{Yeoh}(I_{1})}{\partial
I_{1}},$ (23)

where $I_{1}:={\rm tr}\,{\bf b}$ is the first invariant of the left Cauchy-
Green tensor.

In our model, any deformation transverse to the fibre direction will be
largely due to the base material. Hence we determine the parameters, $c_{k}$,
using experimental data for tissue pulled in the direction transverse to the
fibres. Experimental data is only available for the stress as a function of
stretch in a material; therefore, we use the constitutive relationship (23)
and express the left Cauchy-Green tensor ${\bf b}$ in principal stretches. In
particular, we assume that for uniaxial stretch the principal components of
the stretch, $\lambda_{i}\in\mathbb{R}$, $i=1,\ldots,3$, are given by
$\lambda=\lambda_{1}$ and $\lambda_{2}=\lambda_{3}$, here $\lambda_{1}$ is the
direction of the stretch applied in the experiment. These relations along with
$\lambda_{1}\lambda_{2}\lambda_{3}=1$ (which comes from the incompressibility
assumption) gives us the stress response from the Yeoh model in the principal
spatial direction as

$\displaystyle\sigma_{Yeoh}(\lambda)=2\left(\lambda^{2}-\frac{1}{\lambda}\right)\left[c_{1}+2c_{2}\left(\lambda^{2}+\frac{2}{\lambda}-3\right)+3c_{3}\left(\lambda^{2}+\frac{2}{\lambda}-3\right)\right].$
(24)

The parameters, $c_{k}$, were then obtained via non-linear regression to data
from [47, 19, 3, 31] to determine the muscle material constants and [5] for
the aponeurosis and tendon base material. For the ECM component we utilize
data from a decellularized ECM, which, after decellularization, undergoes
material testing [19]. For the cellular component of the base material,
obtaining mechanical data is more difficult [33]. The cellular component
consists largely of different cells, e.g. muscle fibre cells, satellite cells,
and neuron cell bodies; thus, we can utilize material data obtained from brain
grey matter, which is a mass of cells with little connective tissue [31]. The
fat data is obtained through tensile experiments on subcutaneous fat [3]. This
type of fat is chosen, as it has little anisotropy in the material response,
which is ideal as we model the base material and therefore the fat response as
isotropic. Additionally, this was the only fat data to our knowledge that
contained a tensile response and was in a region near skeletal muscle, as many
other studies investigated the compressive response of fat. The parameter
values from the regression are given in the Supplementary material.

###  4.3 Isochoric response due to along-fibre properties

For the fibre components of the strain-energy function (i.e.
$\Psi_{fibre,tis}$, $tis\in\\{mus,apo,ten\\}$), we use an implicit model where
the constitutive laws are defined in terms of derivatives of these strain-
energies (which is actually what is needed to construct the fictitious Piola-
Kirchhoff stress, $\bar{\bf S}$, see (14)):

$2\bar{I}_{4}\dfrac{\partial\Psi_{fibre,mus}(\bar{I}_{4};\dot{\bar{I}}_{4})}{\partial\bar{I}_{4}}=\sigma_{mus}(\bar{\lambda},\bar{\epsilon}),\quad
2\bar{I}_{4}\dfrac{\partial\Psi_{fibre,tis}(\bar{I}_{4})}{\partial\bar{I}_{4}}=\sigma_{tis}(\bar{\lambda}),$
(25)

for $tis\in\\{apo,ten\\}$, with $\bar{I}_{4}=\bar{\lambda}^{2}$ and
$\dot{\bar{I}}_{4}=2\bar{\lambda}\bar{\epsilon}$ (see subsubsection 2.3.2).
Similar to the formulation of Hill’s model, we can decompose the response from
$\sigma_{mus}$ as

$\sigma_{mus}(\bar{\lambda},\bar{\epsilon})=\sigma_{0}\left(a({\bf
X},t)\,\hat{\sigma}_{act}(\bar{\lambda},\bar{\epsilon})+\hat{\sigma}_{pas}(\bar{\lambda})\right),$
(26)

where $\sigma_{0}>0$ is the maximum isometric stress in the muscle fibres,
$a({\bf X},t)\in[0,1]$ is the fibre activation profile, $\hat{\sigma}_{act}$
is the active stress, and $\hat{\sigma}_{pas}$ is the passive stress in the
fibre. In general, the active stress is defined by two forces of different
nature: elastic forces due to the stretch of the muscle fibres and forces due
to the velocity at which the muscle fibres are being stretched. The active
stress is defined as

$\hat{\sigma}_{act}(\bar{\lambda},\bar{\epsilon}):=\hat{\sigma}_{len}(\bar{\lambda})\,\hat{\sigma}_{vel}(\bar{\epsilon}).$

Here, $\hat{\sigma}_{len}$ is called the force-length function and
$\hat{\sigma}_{vel}$ is the force-velocity function (see Figure 2). For the
sake of clarity, in this case the “velocity” refers to the velocity in which
the fibres are stretched, i.e. the strain rate $\bar{\epsilon}$ (which should
not be confused with either the material or spatial velocity, see (5)). We
also note that the stresses $\hat{\sigma}_{len}$ and $\hat{\sigma}_{vel}$ are
normalized so that $\hat{\sigma}_{len}(1)=1$ and $\hat{\sigma}_{vel}(0)=1$. In
particular, for the quasi-static setting in (18), we have ${\bf v}={\bf 0}$,
and so $\hat{\sigma}_{act}=\hat{\sigma}_{act}(\bar{\lambda})=a({\bf
X},t)\,\hat{\sigma}_{len}(\bar{\lambda})$.

In contrast, aponeurosis and tendon are not subject to active forces, so the
stresses take a simpler form:

$\sigma_{tis}(\bar{\lambda})=\sigma_{0,tis}\hat{\sigma}_{len,tis}(\bar{\lambda}),\quad
tis\in\\{apo,ten\\},$

where $\sigma_{0,tis}>0$ are the maximum isometric stresses and
$\hat{\sigma}_{len,tis}$ are force-length functions of each tissue.

Figure 2: Skeletal muscle material response for the isochoric component of the
base material fit to data from [47] (A). The force-length active and passive
response is fit to data from [92] (B). The response for the force-velocity
component is fit to data from [63] (C).

###  4.4 Parameter fitting for force-length and force-velocity functions

Many different functions have been used to represent the force-length and
velocity relationships [61], including trigonometric (e.g. [75]) and
polynomial functions (e.g. [20]); however, not all these functions are equal,
some will be more computationally efficient than others [61]. A balanced
option developed in previous work [65] uses cubic Bézier curves to capture the
intrinsic properties of muscle. Bézier curves are determined based on a set of
polynomials that smoothly interpolate user defined points [55]. In our case
the user defined points are chosen to achieve a smooth function that
accurately captures the intrinsic muscle properties. A cubic polynomial is
used over higher order polynomials to allow for better local control of the
function.

These Bézier curves have been used in previous one-dimensional modelling
experiments to investigate cyclic work loops [65]; however, in our three-
dimensional model of muscle, they are not as practical. We opted to use a
along fibre response that was twice continuously differentiable. The reason
for this is that that in practice more computational time is required to
compute the resulting stress with the Bézier curve. A cubic polynomial was fit
to the existing Bézier curves from [65] for the passive force-length and
force-velocity response, but with the additional requirement that the function
is twice continuously differentiable everywhere. For the active force-length
relationship a trigonometric function was used opposed to the polynomial fit,
but with the same constraint that the function is twice continuously
differentiable everywhere.

##  5 Numerical experiments

We now present a number computational studies based on our model, exploring
the impact of muscle mass on dynamics, as well as probing the role of
structural changes during cerebral palsy. In subsection 5.2 we demonstrate the
dynamical behavior of the model, contrasting with the use of a quasi-static
approximation. Furthermore, in subsection 5.4 we show an application of the
quasi-static model to muscle effected by cerebral palsy, which demonstrates
the ability of the model to capture the influence of muscular disorders. Our
framework allows us to include patient-specific meshes based on MRI data, and
an experiment demonstrating this is reported in subsection 5.5. All
experiments have been implemented in the finite element library deal.ii
version 8.5.1 [4]. In each, we use mapped hexahedral elements with smallest
$\Delta x=.75$mm, and a time step of $10^{-5}mm$.

###  5.1 Initial and boundary conditions

In what follows, we assume that at time $t=0$, the muscle is at rest,
unactivated, and in an underformed state. Therefore ${\bf u}({\bf
X},0)=\mathbf{0}$, ${\bf v}({\bf X},0)=0$ and ${\bf F}({\bf X},0)={\bf I}$.
The latter implies $D({\bf X},0)=1$ and

$\displaystyle p({\bf
X},0)=\Psi_{vol}^{\prime}(D)\big{|}_{t=0}=\frac{\kappa}{2}\left.\left(D-\frac{1}{D}\right)\right|_{t=0}=0.$

Note that this choice of initial conditions is compatible with the definitions
introduced in section 2.

The boundary conditions that will be implemented in the following numerical
experiments are the same with exception of the exact form of the Dirichlet
condition applied to the muscle. Figure 3 shows diagrams of the geometries in
consideration. We prescribe a zero Dirichlet condition on the $-x$-face of the
geometries (i.e. ${\bf u}_{D}={\bf 0}$). In turn, the condition for the
$+x$-face has the form ${\bf
u}_{D}(\cdot,t)=\left(d(t),0,0\right)^{\texttt{t}}$, with $d(t)$ being a
smooth function that varies according to the experiment. The rest of the
boundary correspond to traction-free surfaces.

$L_{mus}$ 52.0008 [mm] $W_{mus}$ 13.7500 [mm] $H_{mus}$ 5.5783 [mm]
$L_{apo}$ 52.0008 [mm] $\lambda_{0}$ 16.2500 [mm] $\theta_{0}$ 20∘ $L_{mus}$
67.5000 [mm] $T_{apo}$ 0.7500 [mm] $W_{mus}$ 13.7500 [mm]

Figure 3: Geometries considered for the numerical experiments: a block of pure
muscle (top) and the simplified version of the medium gastrocnemius from [64]
(bottom). Note that in the block geometry, the fibres run parallel to the
$x$-axis, whereas in the medium gastrocnemius they run at an angle
$\theta_{0}$ from the $x$-axis. In this case, the relationship between the
length of the muscle $L_{mus}$ and the initial fibre length $\lambda_{0}$ is
given by $\lambda_{0}\sin(\theta_{0})=L_{mus}\sin(\theta_{0}-\gamma_{0})$.

###  5.2 Dynamic muscle mechanics

We begin by discussing the simple experiment discussed in subsubsection 2.3.6.
Consider a block of pure muscle, clamped at one end with fibres oriented along
the x-axis, as in Figure 3. The length of the muscle block is initially
$L_{\mathrm{m}us}$, and points on the $x_{end}$ face of of block are subject
to a Dirichlet boundary condition ${\bf u}_{d}=(d(t),0,0)$, where

$d(t)=\begin{cases}0,&0\leq t\leq 0.05\\\ (t-0.05)\cdot
0.1L_{\mathrm{m}us},&0.05\leq t\leq.015\\\ 0.1L_{\mathrm{m}us}&t\geq
0.15.\end{cases}$

On the transverse faces, a zero traction is prescribed. We simulate the
elastic response in the muscle block over time using a ${\bf
Q}_{2}$-$P_{1}$-$P_{1}$ discretization on a fine mesh.

In Figure 4, we present time traces of the x-components of the displacement.
These are measured at 4 distinct points along the centerline of the muscle
block: at $x(L)$ (the face being pulled), $x_{1}$ close to this face,
$x_{mid}$ at the midpoint of the block and $x_{2}$ close to the clamped face.
The left plot in Figure 4 shows these traces when we use a quasistatic
approximation, and on the right when we use the full dynamic model. We note
that, as expected, a quasistatic approximation leads to all points along the
muscle reacting instantaneously to the pulling (albeit with differing
magnitudes). On the other hand, inertial effects are clearly visible in the
dynamic calculation; an elastic wave propagates towards the clamped end.

Figure 4: The x-component of displacement as a function of time, at different
locations along a parallel-fibered muscle block. One face is fixed, and the
opposite one is initially held at a fixed distance, then pulled with a
constant velocity, and then held fixed again. Right: Quasi-static
calculations. Left: Fully dynamic calculations. Elastic waves are seen clearly
propagating in the dynamic model. The normalization factor is $L_{mus}$ (see
Figure 3).

###  5.3 Isokinetic force-velocity experiments

The force-velocity relationship is a fundamental principle of skeletal muscle
physiology based on Hill’s studies on isolated frog muscles ([2]). In this
study, Hill found that when a muscle is stimulated isometrically and then
suddenly released under a load, the shortening speed increases with decreasing
loads in a hyperbolic relationship ([25]). While Hill’s experiments have been
instrumental to our understanding of the muscle’s force-velocity relationship,
due to the small scale of these experiments, important factors such as mass
and the three-dimensional behavior of muscle tissue were neglected. For
instance, consider single-fibre experiments (or those on small muscles) where
one end is clamped, and the other is moving at a fixed applied velocity of
${\bf v}_{appl}$. One reasonable assumption may be that each point along a
fibre experiences a stretching at a rate given by ${\bf v}_{appl}\cdot{\bf
a}$, where ${\bf a}=\bar{\bf F}{\bf a}_{0}$ is the local direction of the
fibre. We term this an isokinetic assumption; here, mass is included in the
dynamic model. Another reasonable assumption may be that it suffices to
consider a quasi-static approximation (without mass), and in this case the
free end moves with ${\bf u}_{D}=\delta t\,{\bf v}_{appl}$. In a fully dynamic
calculation, the local fibre strain rate is given by (2).

In this experiment, we compare these three models. All computations are
performed on a sample medium gastrocnemius (bottom figure in Figure 3), in
which the fibres are pennate (i.e., not parallel to the direction of the
aponeurosis). The configuration is held fixed along the $-x$-face of the
bottom aponeurosis. The $+x$-face of the top aponeurosis is initially held
fixed while the muscle fibres are activated. Then this face is subject to a
deformation ${\bf u}_{D}=-{\bf v}_{appl}\cdot{\rm e}_{1}$, corresponding to a
shortening.

In each of these simulations, the muscle is activated initially, and then
shortened. In Figure 5, we show the deformed muscles; the colors indicate the
(normalized) z-component of displacement. All three computations yield similar
results during the initial activation phase (i.e., for a low activation
level), but substantially different behaviors once the muscle is shortening
(bottom 3 figures in Figure 5). Once mass effects are included (in the
isokinetic and the dynamic models), the muscle is seen to curve more towards
the end where the shortening is being applied; both the magnitude of this
curvature as well as the pattern of deformation are clearly not the same.

Figure 5: Left to right: A quasi-static, iso-kinetic and fully dynamic
simulation on a model gastrocnemius. Top row: 45% activation. Bottom row:
Activated muscle after shortening phase. The normalization factor corresponds
to the total height of the unit, i.e. $\lambda_{0}\sin(\theta_{0})+2T_{apo}$
(see Figure 3).

###  5.4 Quasi-static muscle mechanics: modelling cerebral palsy

The mechanics of skeletal muscle can be substantially altered by muscular
disorders, aging, and disuse [42, 28]. The altered mechanics come as a result
of changes to the microstructure of the muscle, such as increased
extracellular-matrix (ECM) volume fraction, and changes to the overall muscle
architecture. Understanding how individual changes affect muscle function is
difficult to do using experimental techniques; however, utilizing a modelling
approach we can understand the direct consequence of changes occurring during
muscular disorders. In this work we demonstrate an example of cerebral palsy
(CP) affected muscle. To consider the effects from the micro-mechanical
effects of CP, we vary the volume fraction and stiffness parameters in (9).
Another characteristic of CP affected muscle is increased length of the
sarcomeres; mathematically, this can be modelled as a shift in the fibre
force-length relationship
$\sigma_{mus}(\bar{\lambda}+c_{\text{sarco}},\bar{\epsilon})$ [34].

To evaluate the influence of the CP properties in skeletal muscle, we can
perform a force-length test, using the quasi-static model of muscle (see e.g.
[41]), aiming to mimic typical force-length experiments holding the end of a
muscle fixed. The goal is to study the impact of various structural changes to
a passively-stretched muscle. To construct a passive force length curve, the
muscle is pulled to a given length. The magnitude of the force is then
measured on the end $+x$-face of the muscle. The active force-length curve can
similarly be created by pulling the muscle to a given length, then activating
using a linear activation time ramp to 100%.

We see that our model behaves as expected for typically developed in both the
passive and active setting Figure 6. To investigate the passive mechanics of
CP, we look at the effects of varying $c_{\text{sarco}}$ and $\alpha$
independently. We note this cannot be done in vivo, as both of these
properties vary simultaneously during CP [79]. Using our model we are able to
determine that there is a larger effect from the ECM relative to the sarcomere
length in the passive setting (Figure 6). The benefit of our three dimensional
continuum model is that we can customize the sample geometry to more
accurately capture the deformation of muscle. In Figure 6 (B), we utilize a
pennate muscle geometry (Figure 6 D), where the muscle fibres are orientated
on an angle and connected to a thin stiff sheet of aponeurosis material. Using
this geometry allows us to investigate to common architectural changes in CP
muscle, such the decrease in physiological cross-sectional area and volume
[6], on the muscle mechanics. We ask the question: which of the changes –
those to the physiological cross-section, or to the properties of the skeletal
muscle itself – impact the force output more? To answer this, we independently
varied the physiological cross section and the tissue properties. In the
computational experiments corresponding to typical muscle (TD) there is an
$\approx 30\%$ decrease in force when the physiological cross section is
decreased by 30%.(Figure 6). This is to be expected, since it is often assumed
that physiological cross section is proportional to muscle force [39]. When
looking at the CP muscle tissue we instead find that there is less effect from
the physiological cross section ($\approx 22\%$) (Figure 6), which indicates
an influence from the CP material properties. Our model is able to capture the
influence of changes to the material properties and architecture of CP muscle,
and determine their influence on the overall muscle mechanics. Due to the
highly complex structure and variability that occurs in CP affected muscle,
this can not be done using experimental techniques.

Figure 6: Modelling the influence of CP on muscle mechanics. The passive
force-stretch relationship is shown in (A) for a block of muscle tissue (C).
The sarcomere length was shifted by a factor a $c_{\text{sarco}}$ and the
volume fraction of the ECM was varied to represent the effects of CP. Dotted
lines represent the curves corresponding to a CP base material. The active
force-stretch relationship is shown in (B) for a pennate muscle geometry (D).
Here, we show the influence of CP muscle geometries on the total force
produced by muscle. The dotted lines correspond to CP muscle material
properties ($\alpha=0.4$, $\beta=0.2$, and $c_{\text{sarco}}=0.0$), while the
solid lines correspond to TD muscle properties ($\alpha=0.02$, $\beta=0.1$,
and $c_{\text{sarco}}=0.0$)

###  5.5 MRI-derived geometry

We can utilize a realistic MRI-derived geometry with our model to solve our
model for subject-specific data (Figure 7), which were obtained by D’Souza et
al. [16]. Using these data allows us to account for variations in the
structure of muscle that may vary between subjects. This could include changes
to muscle thickness, length, and volume, all of which will have an effect on
the mechanics of skeletal muscle [40]. Figure 7 shows a simulation where the
muscle was stretched in the $x$ direction to a whole muscle stretch of 1.1,
and subsequently activated to 100% activation at a fixed length. The $+z$ and
$-z$ faces have a portion of the surface covered by an aponeurosis, which the
muscle fibres insert into. The fibres in this model run at 25 degrees from the
$x$ axis.

Figure 7: Magnitude of displacement (Top) and $z$ direction displacement at
100% activation on a subject specific MRI-derived geometry. Displacements
shown are in units of meters.

In order to use our computational framework with this more realistic geometry
we need the information from post-processed MRI scans to be able to define the
mesh of such geometries. We also need DTI scans to be able to map the muscle
fibre orientation within such a geometry; for more details see, for instance
[88, 12]. As one can imagine, the fibre directions inside a muscle can be
quite complex. Fortunately, we only need a snapshot of the geometry and the
fibre orientation vectors to be used at the beginning of the simulation, as
part of the initial setup of the system.

###  5.6 Validation of models

A common (and legitimate) criticism expressed by scientists in physiology of
complex mathematical models concerns their validity: if many parameters are
involved in the model which are fit to experimental data, how do we guard
against over-fitting? In our framework, we have fit parameters describing
tissue properties. We can then compare quantities such as displacement, fibre
orientations, energies, and forces which are not fit; these comparisons are
with in situ experiments. Naturally, our simulations cannot reproduce the
exact values seen in experiment since we necessarily work with idealized
geometries and configurations. We can, however, compare qualitative trends
between experimental data and computational results.

As a first such validation, in [88] we compared MRI data on muscle bulging
during fixed-end contractions, with a simulation based on our model on an MRI-
derived mesh. Notably, even though the location of the aponeurosis in our
computation was only indirectly inferred from imaging data, and though our
computations did not include the effect of tendons or adjacent muscles, our
simulations predicted magnitudes of both inward and outward bulging which were
consistent with those experimentally measured (Figure 12 in [88].The data we
used to obtain experimental parameters did not include muscle bulging
information.

We also found similar patterns of anisotropy in transverse bulging as previous
ultrasound measures [60]. Further, an ultrasound study of human ankle
plantarflexors showed that muscles with lower initial fibre pennation angles
(less than 15 degrees) tend to bulge, whereas more pennate muscles tend to
thin as they activate [59], a feature of muscle contraction that we were able
to replicate with our model in [88].

Next, in [71], we considered fixed-end contractions in a parallel-fibered
muscle which was subject to a compressive force transverse to the fibers, and
the fibres are then activated. Our results showed that before activation,
compression would cause muscle thickness to decrease and the fibres to change
their pennation angle; this matched the trends reported in [72]. Our
simulations also revealed a complex dependance of the longtitudinal force
(oriented between the fixed ends), the length of the muscle block, and the
initial fibre pennation angle. The trend we observed: compressive transverse
loading lead to an increased longtitudinal force at longer lengths and a
reduction at shorter lengths, consistent with experiments reported by [78].

Finally, to investigate the importance of mass and muscle size on the
dynamics, in situ experiments were conducted on the rat plantaris muscle, in
which the muscle was subjected to cyclic contractions, coupled with neuronal
excitations to cause activation in synch [67]. Accelerations in the $x$
direction, $a_{mid}$ and $a_{end}$ were measured at the midpoint $L_{mus}/2$
and at the end point of the muscle $L_{mus}$, respectively. If the muscle
mechanics of the rat plantaris muscle were genuinely well-approximated by a
single fibre quasi-static Hill-type model, the scaled quantities
$2a_{mid}/(L_{mus})$ and $a_{end}/L_{mus}$ should be the same. Experimentally,
these are not; the difference is more pronounced when mass was added to the
muscle, and also in larger muscles.

We ran simulations with the fully dynamic model on a muscle-aponeurosis
configuration as in Figure 3, and subjected the $+x$-face of the top
aponeurosis to cyclic contractions; simultaneously, the muscle fibres were
activated in a synchronized manner. We found comparable patterns of tissue
accelerations in the longitudinal direction across the dynamic model during
cyclic contractions [64] to that of the experiments on rat plantaris muscle
[67], again scaling the accelerations by a length scale corresponding to the
location in the muscle. We found a larger difference in scaled acceleration in
the middle of the muscle relative to the end with greater muscle size for the
model and greater added mass for the in situ muscle.

##  6 Discussion and validation

In this paper we proposed a mathematical model to describe the three
dimensional deformations of skeletal muscles. We used the principles of
continuum mechanics to derive this model. Muscle is considered as a nearly
incompressible, fibre-reinforced and transversely isotropic composite
material; internal chemical reactions are triggered to produce forces which
deform the tissues without any external work done on the tissues, while
passive elastic properties in the muscle fibres resist the deformation in the
along fibre direction. The idea of a base material was used to encapsulate the
contribution from several other tissues and fluids surrounding the fibres:
blood, water, extra- and intra-cellular matrix forms part of this material.
Experimental data was used to fit parameters in the model.

We also proposed a numerical approach to approximate the deformations of the
tissues. A semi-implicit scheme is used to discretize the time variables
whereas a conforming finite element method is used to discretize the spatial
variables. As a result of such time stepping, a CFL condition governs both the
time step size and mesh size in order for the computation to be stable. Note
that the implementation of our method was not optimized to match experimental
data. Instead, the outputs of our model are a result of the assumptions of the
physics and mechanics in our model and its implementation.

Finally, we remark that many natural processes are not considered in this
study. For example, muscle contraction in living beings produces heat. Also,
fluids as water and blood form part of the overall composition of muscle
tissues. Even though the influence of fluids is implicitly added to our model
through the experimental data we used to fit the model parameters, the current
model does not explicitly consider the effects of internal fluids as part of
the deformation of the tissues. Finally, chemical reactions leading to the
activation of muscle fibres is not included explicitly; only the action of
these reactions is encoded in our model through the active stress
$\hat{\sigma}_{act}$ and the form of the activation function which in turn
depends on neuronal excitation, and can be modeled via a system of ODEs. These
features correspond to an area of much current research which are not included
in our model, and are left as future work.

## Acknowledgments

Sebastián A. Domínguez-Rivera thanks the financial support of PIMS-Canada
through a PIMS postdoctoral fellowship. Nilima Nigam thanks the support of
NSERC-Canada through the Discovery program.

##### Note for the reader

Please note that this manuscript represents a preprint only and has not been
(or is in the process of being) peer-reviewed. A DOI link will be made
available for this ArXiv preprint as soon as the peer-reviewed version is
published online.

## References

  * [1] T. Ahamed, M. B. Rubin, B. A. Trimmer, and L. Dorfmann.  Time-dependent behavior of passive skeletal muscle.  Continuum Mechanics and Thermodynamics, 28(1-2):561–577, 2016.
  * [2] J. Alcazar, R. Csapo, I. Ara, and L. M. Alegre.  On the shape of the force-velocity relationship in skeletal muscles: The linear, the hyperbolic, and the double-hyperbolic.  Front. Physiol. 10:769, 2019.
  * [3] N. Alkhouli, J. Mansfield, E. Green, J. Bel, B. Knight, N. Liversedge, J. C. Tham, R. Welbourn, A. C. Shore, K. Kos, and C. P. Winlove.  The mechanical properties of human adipose tissues and their relationships to the structure and composition of the extracellular matrix.  American Journal of Physiology - Endocrinology and Metabolism, 305(12):1427–1435, 2013.
  * [4] D. Arndt, W. Bangerth, D. Davydov, T. Heister, L. Heltai, M. Kronbichler, M. Maier, J.-P. Pelteret, B. Turcksin, and D. Wells.  The deal.II library, version 8.5.  Journal of Numerical Mathematics, 25(3):137–146, 2017.
  * [5] E. Azizi, G. M. Halenda, and T. J. Roberts.  Mechanical properties of the gastrocnemius aponeurosis in wild turkeys.  Integrative and Comparative Biology, 49(1):51–58, 2009.
  * [6] R. S. Barrett and G. A. Lichtwark.  Gross muscle morphology and structure in spastic cerebral palsy: a systematic review.  Developmental Medicine and Child Neurology, 52(9):794–804, 2010\.
  * [7] R. J. Baskin and P. J. Paolini.  Volume change and pressure development in muscle during contraction.  American Journal of Physiology-Legacy Content, 213(4):1025–1030, 1967.
  * [8] C. Bleiler, P. Ponte Castañeda, and O. Röhrle.  A microstructurally-based, multi-scale, continuum-mechanical model for the passive behaviour of skeletal muscle tissue.  Journal of the Mechanical Behavior of Biomedical Materials, 97(April):171–186, 2019.
  * [9] S. S. Blemker and S. L. Delp.  Three-dimensional representation of complex muscle architectures and geometries.  Annals of Biomedical Engineering, 33(5):661–673, 2005.
  * [10] M. Böl and S. Reese.  Micromechanical modelling of skeletal muscles based on the finite element method.  Computer Methods in Biomechanics and Biomedical Engineering, 11(5):489–504, 2008.
  * [11] L. M. Bollinger.  Potential contributions of skeletal muscle contractile dysfunction to altered biomechanics in obesity.  Gait & Posture, 56:100–107, 2017.
  * [12] B. Bolsterlee.  A new framework for analysis of three-dimensional shape and architecture of human skeletal muscles from in vivo imaging data.  Journal of Applied Physiology, 132(3):712–725, 2022.  PMID: 35050794.
  * [13] T. J. Burkholder and R. L. Lieber.  Sarcomere length operating range of vertebrate muscles during movement.  J. Exp. Biol., 204(9):1529–1536, 2001.
  * [14] S. L. Delp, F. C. Anderson, A. S. Arnold, P. Loan, A. Habib, C. T. John, E. Guendelman, and D. G. Thelen.  OpenSim: Open-source software to create and analyze dynamic simulations of movement.  IEEE Transactions on Biomedical Engineering, 54(11):1940–1950, 2007\.
  * [15] T. J. M. Dick, A. A. Biewener, and J. M. Wakeling.  Comparison of human gastrocnemius forces predicted by Hill-type muscle models and estimated from ultrasound images.  Journal of Experimental Biology, 220(9):1643–1653, 2017.
  * [16] A. D’Souza, B. Bolsterlee, A. Lancaster, and R. D. Herbert.  Muscle architecture in children with cerebral palsy and ankle contractures: an investigation using diffusion tensor imaging.  Clinical Biomechanics, 68(June):205–211, 2019.
  * [17] J. C. Gardiner and J. A. Weiss.  Simple shear testing of parallel-fibered planar soft tissues.  Journal of Biomechanical Engineering, 123(2):170–175, 2001.
  * [18] M. H. Gfrerer and B. Simeon.  Fiber-based modeling and simulation of skeletal muscles.  Multibody System Dynamics, 52(1):1–30, 2021.
  * [19] A. R. Gillies, L. R. Smith, R. L. Lieber, and S. Varghese.  Method for decellularizing skeletal muscle without detergents or proteolytic enzymes.  Tissue Engineering - Part C: Methods, 17(4):383–389, 2011.
  * [20] J. He, W. S. Levine, and G. E. Loeb.  Feedback Gains for Correcting Small Perturbations to Standing Posture.  IEEE Transactions on Automatic Control, 36(3):322–332, 1991.
  * [21] S. Hedenstierna, P. Halldin, and K. Brolin.  Evaluation of a combination of continuum and truss finite elements in a model of passive and active muscle tissue.  Computer Methods in Biomechanics and Biomedical Engineering, 11(6):627–639, 2008.  PMID: 18642161.
  * [22] T. Heidlauf, T. Klotz, C. Rode, T. Siebert, and O. Röhrle.  A continuum-mechanical skeletal muscle model including actin-titin interaction predicts stable contractions on the descending limb of the force-length relation.  PLOS Computational Biology, 13(10):1–25, 10 2017.
  * [23] T. Heidlauf and O. Röhrle.  A multiscale chemo-electro-mechanical skeletal muscle model to analyze muscle contraction and force generation for different muscle fiber arrangements.  Frontiers in Physiology, 5, 2014.
  * [24] W. Herzog.  The multiple roles of titin in muscle contraction and force production.  Biophysical Reviews, 10(4):1187–1199, 2018.
  * [25] A. V. Hill.  The heat of shortening and the dynamic constants of muscle.  Proceedings of the Royal Society of London. Series B-Biological Sciences, 126(843):136–195, 1938.
  * [26] N. C. Holt, J. M. Wakeling, and A. A. Biewener.  The effect of fast and slow motor unit activation on whole-muscle mechanical performance: the size principle may not pose a mechanical paradox.  Proceedings of the Royal Society B: Biological Sciences, 281(1783):20140002, 2014.
  * [27] G. A. Holzapfel.  Nonlinear solid mechanics: a continuum approach for engineering, 2nd edn chichester.  UK: John Wiley & Sons, 2001.
  * [28] J. J. Howard and W. Herzog.  Skeletal Muscle in Cerebral Palsy: From Belly to Myofibril.  Frontiers in Neurology, 12(February):1–15, 2021.
  * [29] Huxley and A. F.  Muscle structure and theories of contraction.  Prog. Biophys. Mol. Biol., 7:255–318, jan 1957.
  * [30] H. E. Huxley.  Past, present and future experiments on muscle.  Philosophical Transactions of the Royal Society of London. Series B: Biological Sciences, 355(1396):539–543, apr 2000.
  * [31] X. Jin, F. Zhu, H. Mao, M. Shen, and K. H. Yang.  A comprehensive experimental study on material properties of human brain tissue.  Journal of Biomechanics, 46(16):2795–2801, 2013.
  * [32] K. R. Knaus, G. G. Handsfield, and S. S. Blemker.  A 3d model of the soleus reveals effects of aponeuroses morphology and material properties on complex muscle fascicle behavior.  Journal of Biomechanics, 130:110877, 2022.
  * [33] R. N. Konno, N. Nigam, and J. M. Wakeling.  Modelling extracellular matrix and cellular contributions to whole muscle mechanics.  PLoS ONE, 16(4 April 2021):1–20, 2021.
  * [34] R. N. Konno, N. Nigam, J. M. Wakeling, and S. A. Ross.  The Contributions of Extracellular Matrix and Sarcomere Properties to Passive Muscle Stiffness in Cerebral Palsy.  Frontiers in Physiology, 12(January):1–13, 2022.
  * [35] M. N. Korta, D. Nakhostin, L. Ruby, T. Frauenfelder, M. Rominger, and S. Sanabria.  Speed of sound and shear wave speed for calf soft tissue composition and nonlinearity assessment.  Quant Imaging Med Surg., 11(9), 2021.
  * [36] J. Lamsfuss and S. Bargmann.  Skeletal muscle: Modeling the mechanical behavior by taking the hierarchical microstructure into account.  Journal of the Mechanical Behavior of Biomedical Materials, 122:104670, 2021.
  * [37] J. Li.  Development and validation of a finite-element musculoskeletal model incorporating a deformable contact model of the hip joint during gait.  Journal of the Mechanical Behavior of Biomedical Materials, 113:104136, 2021.
  * [38] J. Li, Y. Lu, S. C. Miller, Z. Jin, and X. Hua.  Development of a finite element musculoskeletal model with the ability to predict contractions of three-dimensional muscles.  Journal of Biomechanics, 94:230–234, 2019.
  * [39] R. L. Lieber and J. Fridén.  Functional and clinical significance of skeletal muscle architecture.  Muscle & Nerve: Official Journal of the American Association of Electrodiagnostic Medicine, 23(11):1647–1666, 2000.
  * [40] R. L. Lieber and J. Fridén.  Functional and clinical significance of skeletal muscle architecture.  Muscle Nerve, 23(11):1647–1666, 2000.
  * [41] R. L. Lieber, S. Steinman, I. A. Barash, and H. Chambers.  Structural and functional changes in spastic skeletal muscle, may 2004\.
  * [42] R. L. Lieber and S. R. Ward.  Cellular mechanisms of tissue fibrosis. 4. structural and functional consequences of skeletal muscle fibrosis.  American Journal of Physiology-Cell Physiology, 305(3):C241–C252, 2013.
  * [43] Y. T. Lu, L. Beldie, B. Walker, S. Richmond, and J. Middleton.  Parametric study of a hill-type hyperelastic skeletal muscle model.  Proceedings of the Institution of Mechanical Engineers, Part H: Journal of Engineering in Medicine, 225(5):437–447, 2011.  PMID: 21755774.
  * [44] S. P. Magnusson, P. Hansen, P. Aagaard, J. Brønd, P. Dyhre-Poulsen, J. Bojsen-Moller, and M. Kjaer.  Differential strain patterns of the human gastrocnemius aponeurosis and free tendon, in vivo.  Acta Physiologica Scandinavica, 177(2):185–195, 2003.
  * [45] B. Maier and M. Schulte.  Mesh generation and multi-scale simulation of a contracting muscle–tendon complex.  Journal of Computational Science, 59:101559, 2022.
  * [46] A. L. Mescher.  Muscle Tissue.  McGraw-Hill Education, New York, NY, 2018.
  * [47] M. Mohammadkhah, P. Murphy, and C. K. Simms.  The in vitro passive elastic response of chicken pectoralis muscle to applied tensile and compressive deformation.  Journal of the Mechanical Behavior of Biomedical Materials, 62:468–480, 2016.
  * [48] E. K. Moo, T. R. Leonard, and W. Herzog.  The sarcomere force–length relationship in an intact muscle–tendon unit.  Journal of Experimental Biology, 223(6), 03 2020.  jeb215020.
  * [49] N. Nigam, J. M. Wakeling, and H. Rahemi.  Regionalization muscle activity causes changes to the magnitude and direction of the force from the whole muscles– a modelling study.  Frontiers in physiology, 5:298, 2014.
  * [50] R. W. Ogden.  Non-Linear Elastic Deformations.  Dover Civil and Mechanical Engineering. Dover Publications, 1997.
  * [51] C. W. J. Oomens, M. Brekelmans, and F. P. Baaijens.  Biomechanics: Concepts and Computation.  Cambridge Texts in Biomedical Engineering. Cambridge University Press, 2009.
  * [52] C. W. J. Oomens, M. Maenhout, C. H. Van Oijen, M. R. Drost, and F. P. Baaijens.  Finite element modelling of contracting skeletal muscle.  Philosophical Transactions of the Royal Society of London. Series B: Biological Sciences, 358(1437):1453–1460, 2003.
  * [53] J.-P. Pelteret and A. McBride.  The deal.ii tutorial step-44: Three-field formulation for non-linear solid mechanics, 2012.
  * [54] E. J. Perreault, C. J. Heckman, and T. G. Sandercock.  Hill muscle model errors during movement are greatest within the physiologically relevant range of motor unit firing rates.  Journal of Biomechanics, 36(2):211–218, 2003.
  * [55] H. Prautzsch.  Bézier and B-spline techniques / Hartmut Prautzsch, Wolfgang Boehm, Marco Paluszny. Mathematics and visualization. Springer, 2002.
  * [56] P. P. Purslow and J. A. Trotter.  The morphology and mechanical properties of endomysium in series-fibred muscles: variations with muscle length.  Journal of Muscle Research & Cell Motility, 15(3):299–308, 1994\.
  * [57] H. Rahemi, N. Nigam, and J. M. Wakeling.  Regionalizing muscle activity causes changes to the magnitude and direction of the force from whole muscles-a modeling study.  Frontiers in Physiology, 5 AUG(August):1–10, 2014.
  * [58] H. Rahemi, N. Nigam, and J. M. Wakeling.  The effect of intramuscular fat on skeletal muscle mechanics: implications for the elderly and obese.  Journal of The Royal Society Interface, 12(109):20150365, 2015.
  * [59] A. Randhawa, M. E. Jackman, and J. M. Wakeling.  Muscle gearing during isotonic and isokinetic movements in the ankle plantarflexors.  European Journal of Applied Physiology, 113(2):437–447, 2013.
  * [60] A. Randhawa and J. M. Wakeling.  Transverse anisotropy in the deformation of the muscle during dynamic contractions.  Journal of Experimental Biology, 221(15):jeb175794, 2018.
  * [61] R. Rockenfeller and M. Günther.  How to model a muscle’s active force–length relation: A comparative study.  Computer Methods in Applied Mechanics and Engineering, 313:321–336, jan 2017.
  * [62] O. Röhrle and A. J. Pullan.  Three-dimensional finite element modelling of muscle forces during mastication.  Journal of Biomechanics, 40(15):3363–3372, 2007.
  * [63] H. Roots, G. W. Offer, and K. W. Ranatunga.  Comparison of the tension responses to ramp shortening and lengthening in intact mammalian muscle fibres: Crossbridge and non-crossbridge contributions.  Journal of Muscle Research and Cell Motility, 28(2-3):123–139, feb 2007.
  * [64] S. A. Ross, S. Domínguez, N. Nigam, and J. M. Wakeling.  The Energy of Muscle Contraction. III. Kinetic Energy During Cyclic Contractions.  Frontiers in Physiology, 12(April):1–16, 2021.
  * [65] S. A. Ross, N. Nigam, and J. M. Wakeling.  A modelling approach for exploring muscle dynamics during cyclic contractions.  PLoS computational biology, 14(4):e1006123, 2018.
  * [66] S. A. Ross, N. Nigam, and J. M. Wakeling.  A modelling approach for exploring muscle dynamics during cyclic contractions.  PLoS Computational Biology, 14(4):e1006123, 2018.
  * [67] S. A. Ross, B. Rimkus, N. Konow, A. A. Biewener, and J. M. Wakeling.  Added mass in rat plantaris muscle causes a reduction in mechanical work.  Journal of Experimental Biology, 223(19), 10 2020.  jeb224410.
  * [68] S. A. Ross, D. S. Ryan, S. Domínguez, N. Nigam, and J. M. Wakeling.  Size, history-dependent, activation and three-dimensional effects on the work and power produced during cyclic muscle contractions.  Integrative and Comparative Biology, 58(2):232–250, 2018.
  * [69] S. A. Ross and J. M. Wakeling.  Muscle shortening velocity depends on tissue inertia and level of activation during submaximal contractions.  Biology letters, 12(6):20151041, 2016.
  * [70] A. Roux, J. Lecompte, I. Iordanoff, and S. Laporte.  Modeling of muscular activation of the muscle-tendon complex using discrete element method.  Computer Methods in Biomechanics and Biomedical Engineering, 24(11):1184–1194, 2021.  PMID: 33416406.
  * [71] D. S. Ryan, S. Domínguez, S. A. Ross, N. Nigam, and J. M. Wakeling.  The energy of muscle contraction. ii. transverse compression and work.  Frontiers in Physiology, 11:1392, 2020.
  * [72] D. S. Ryan, N. Stutzig, T. Siebert, and J. M. Wakeling.  Passive and dynamic muscle architecture during transverse loading for gastrocnemius medialis in man.  Journal of Biomechanics, 86:160–166, 2019.
  * [73] O. Röhrle, J. Davidson, and A. Pullan.  A physiologically based, multi-scale model of skeletal muscle structure and function.  Frontiers in Physiology, 3, 2012.
  * [74] S. Schmitt, M. Günther, and D. F. B. Häufle.  The dynamics of the skeletal muscle: A systems biophysics perspective on muscle modeling with the focus on hill-type muscle models.  GAMM-Mitteilungen, 42(3):e201900013, 2019.
  * [75] S. H. Scott and D. A. Winter.  A comparison of three muscle pennation assumptions and their effect on isometric and isotonic force.  Journal of Biomechanics, 24(2):163–167, jan 1991.
  * [76] B. Sharafi and S. S. Blemker.  A micromechanical model of skeletal muscle to explore the effects of fiber and fascicle geometry.  Journal of Biomechanics, 43(16):3207–3213, 2010.
  * [77] J. C. Simo, R. L. Taylor, and K. S. Pister.  Variational and projection methods for the volume constraint in finite deformation elasto-plasticity.  Computer Methods in Applied Mechanics and Engineering, 51(1):177 – 208, 1985.
  * [78] D. A. Sleboda and T. J. Roberts.  Internal fluid pressure influences muscle contractile force.  Proceedings of the National Academy of Sciences, 117(3):1772–1778, 2020.
  * [79] L. R. Smith, L. Gerace-Fowler, and R. L. Lieber.  Muscle extracellular matrix applies a transverse stress on fibers with axial strain.  Journal of Biomechanics, 44(8):1618–1620, 2011.
  * [80] L. A. Spyrou, M. Agoras, and K. Danas.  A homogenization model of the Voigt type for skeletal muscle.  Journal of Theoretical Biology, 414:50–61, 2017.
  * [81] M. Stefanati, C. Villa, Y. Torrente, and J. F. Rodriguez Matas.  A mathematical model of healthy and dystrophic skeletal muscle biomechanics.  Journal of the Mechanics and Physics of Solids, 134, 2020.
  * [82] C. Tanner, J. A. Schnabel, D. L. Hill, D. J. Hawkes, M. O. Leach, and D. R. Hose.  Factors influencing the accuracy of biomechanical breast models.  Medical physics, 33(6Part1):1758–1769, 2006.
  * [83] M. M. Van Der Krogt, L. Bar-On, T. Kindt, K. Desloovere, and J. Harlaar.  Neuro-musculoskeletal simulation of instrumented contracture and spasticity assessment in children with cerebral palsy.  Journal of NeuroEngineering and Rehabilitation, 13(1):1–11, 2016\.
  * [84] M. Van Loocke, C. G. Lyons, and C. K. Simms.  Viscoelastic properties of passive skeletal muscle in compression: stress-relaxation behaviour and constitutive modelling.  Journal of Biomechanics, 41(7):1555–1566, 2008.
  * [85] K. M. Virgilio, K. S. Martin, S. M. Peirce, and S. S. Blemker.  Multiscale models of skeletal muscle reveal the complex effects of muscular dystrophy on tissue mechanics and damage susceptibility.  Interface Focus, 5(2), 2015.
  * [86] K. Volokh.  Mechanics of soft materials.  Springer, 2016.
  * [87] J. M. Wakeling, S. S. Lee, A. S. Arnold, M. de Boef Miara, and A. A. Biewener.  A muscle’s force depends on the recruitment patterns of its fibers.  Annals of Biomedical Engineering, 40(8):1708–1720, 2012.
  * [88] J. M. Wakeling, S. A. Ross, D. S. Ryan, B. Bolsterlee, R. Konno, S. Domínguez, and N. Nigam.  The energy of muscle contraction. i. tissue force and deformation during fixed-end contractions.  Frontiers in Physiology, 11:813, 2020.
  * [89] J. D. Webb, S. S. Blemker, and S. L. Delp.  3d finite element models of shoulder muscles for computing lines of actions and moment arms.  Computer Methods in Biomechanics and Biomedical Engineering, 17(8):829–837, 2014.  PMID: 22994141.
  * [90] J. Weickenmeier, M. Itskov, E. Mazza, and M. Jabareen.  A physically motivated constitutive model for 3d numerical simulation of skeletal muscles.  International Journal for Numerical Methods in Biomedical Engineering, 30(5):545–562, 2014.
  * [91] J. A. Weiss, B. N. Maker, and S. Govindjee.  Finite element implementation of incompressible, transversely isotropic hyperelasticity.  Computer methods in applied mechanics and engineering, 135(1-2):107–128, 1996.
  * [92] T. M. Winters, M. Takahashi, R. L. Lieber, and S. R. Ward.  Whole muscle length-tension relationships are accurately modeled as scaled sarcomeres in rabbit hindlimb muscles.  Journal of Biomechanics, 44(1):109–115, 2011.
  * [93] P. Wriggers.  Nonlinear Finite Element Methods.  Springer Berlin Heidelberg, 2009.
  * [94] O. H. Yeoh.  Some forms of the strain energy function for rubber.  Rubber Chemistry and Technology, 66(5):754–771, 1993.
  * [95] C. A. Yucesoy, B. H. Koopman, P. A. Huijing, and H. J. Grootenboer.  Three-dimensional finite element modeling of skeletal muscle using a two-domain approach: linked fiber-matrix mesh model.  Journal of Biomechanics, 35(9):1253–62, 2002.
  * [96] F. E. Zajac.  Muscle and tendon: properties, models, scaling, and application to biomechanics and motor control.  Critical Reviews in Biomedical Engineering, 17(4):359–411, 1989\.
  * [97] V. M. Zatsiorsky and B. I. Prilutsky.  Biomechanics of Skeletal Muscles.  Human Kinetics, 2012.
  * [98] J. Östh, K. Brolin, and R. Happee.  Active muscle response using feedback control of a finite element human arm model.  Computer Methods in Biomechanics and Biomedical Engineering, 15(4):347–361, 2012.  PMID: 21294008.

##  7 Supplemental information

In this section we provide the parameters obtained through fitting described
in section 4. For completeness, we also record parameters already reported in
the paper.

###  7.1 Bulk moduli

As described in subsection 4.1, each of the muscle, aponeurosis, tendon and
fat components of an MTU are nearly incompressible. Aponeurosis and tendon are
assumed to have the same bulk modulus, $\kappa_{apo}=\kappa_{ten}=1\times
10^{8}$Pa, [49]. The bulk modulus for extracellular material is
$\kappa_{ecm}=1\times 10^{6}$Pa [9, 57], while that for cellular material
(omitting ECM) and fat are the same $\kappa_{cell}=\kappa_{fat}=1\times
10^{7}$Pa [33]. The bulk modulus of muscle is modeled as a linear combination
of that of fat, cellular material and ECM, depending on their respective
volume fractions.

###  7.2 Yeoh models of elasticity

Each of muscle, tendon and aponeurosis are modelled as fibre-reinforced
composites. The base material is isotropic and nonlinear, and the isochoric
responses as modelled as in section 4. The isochoric base materials for
muscle, aponeurosis and tendon are described using Yeoh-type models; Neo-
hookean models have also been tried and yield similar results [49]. The
parameters $c_{k}$ are obtained via non-linear regression to data from [47,
19, 3, 31] to determine the muscle material constants and [5] for the
aponeurosis and tendon base material, Table 1.

Parameter | Value/Range of Values
---|---
$c_{\text{1,mus}}$ | 3703
$c_{\text{2,mus}}$ | -707.7
$c_{\text{3,mus}}$ | 123.2
$c_{\text{1,apo}}$ | 4.6896264
$c_{\text{2,apo}}$ | -3.455141
$c_{\text{3,apo}}$ | 484.92055
$c_{\text{1,ten}}$ | 4.6896264
$c_{\text{2,ten}}$ | -3.455141
$c_{\text{3,ten}}$ | 484.92055
Table 1: List of the material parameters used in the Yeoh model for the base
material for each of the tissues used in the model.

Based on fits to experimental measurements on human breast tissue [82], the
strain energy of fat is modelled as $\Psi_{fat}(\bar{I_{1}})=0.13\times
10^{6}(\bar{I}_{1}-3)$ Pa, [58].

###  7.3 Along-fibre properties

As discussed in section 4, we initially fit passive, active force-length and
active force-velocity along-fibre properties of muscle to experimental data
using Bezier splines, [66]. These have directly been used in 1-D dynamic
simulations. For the fully 3-dimensional calculations, these expressions are
less convenient. The along-fibre passive response of muscle is modelled via a
piecewise cubic polynomial which is a least-squares best fit to the Bezier
curve and which has smooth derivatives:

$\hat{\sigma}_{\text{pass}}(\lambda)=\left\\{\begin{aligned}
&0&&0\leq\lambda\leq 1.0\\\ &2.353(\lambda-1.0)^{2}&&1.0\leq\lambda\leq
1.25\\\ &3.44(\lambda-1.25)^{2}+1.18(\lambda-1.25)+0.147&&1.25\leq\lambda\leq
1.5\\\ &0.427(\lambda-1.5)^{2}+2.90(\lambda-1.5)+0.656&&1.5\leq\lambda\leq
1.65\\\ &3.023(\lambda-1.65)+1.1&&\lambda>1.65,\end{aligned}\right.$ (27)

The active force-length response for muscle fibres is given by piecewise
trigonometric fit to the Bezier splines in [66].

$\hat{\sigma}_{\text{act}}(\lambda)=\left\\{\begin{aligned}
&(0.642\sin(1.29\lambda+0.629)\\\ &+0.325\sin(5.31\lambda-4.52)\\\
&+0.328\sin(6.74\lambda+1.69)\\\ &+0.015\sin(19.8\lambda-7.39)&&\text{if
}0.4\leq\lambda\leq 1.75\\\ &+0.139\sin(8.04\lambda+2.54)\\\
&+0.0018\sin(32.2\lambda-6.45)\\\ &+0.012\sin(23.2\lambda-2.64))\\\ &&&\\\
0&&&\text{otherwise}.\end{aligned}\right.$ (28)

Both the passive response and the active force-length response are based on
fits to experimental data on rabbit hindlimb muscles, [92].

The active force-velocity relationship for muscle fibres is given by

$\hat{\sigma}_{\text{vel}}(\bar{\epsilon})=\left\\{\begin{aligned}
&0&&0\leq\bar{\epsilon}\leq x_{1}\\\
&0.2579(\bar{\epsilon}-x_{1})^{3}+0.1431(\bar{\epsilon}-x_{1})^{2}&&x_{1}\leq\bar{\epsilon}\leq
x_{2}\\\
&29.8255(\bar{\epsilon}-x_{2})^{3}-0.9435(\bar{\epsilon}-x_{2})^{2}+0.9703(\bar{\epsilon}-x_{2})^{1}+0.3503&&x_{2}\leq\bar{\epsilon}\leq
x_{3}\\\
&-3165.6847(\bar{\epsilon}-x_{3})^{3}+186.1961(\bar{\epsilon}-x_{3})^{2}+6.0908(\bar{\epsilon}-x_{3})^{1}+1&&x_{3}\leq\bar{\epsilon}\leq
x_{4}\\\
&0.6882(\bar{\epsilon}-x_{4})^{3}-1.4139(\bar{\epsilon}-x_{4})^{2}+0.9678(\bar{\epsilon}-x_{4})^{1}+1.3743&&x_{4}\leq\bar{\epsilon}\leq
x_{5}\\\ &1.5950&&\bar{\epsilon}\geq x_{5},\end{aligned}\right.$ (29)

with $x_{1}=-1.2\bar{\epsilon}_{0}$,
$x_{2}=-0.25\bar{\epsilon}_{0}$,$x_{3}=0$, $x_{4}=0.05\bar{\epsilon}_{0}$, and
$x_{5}=0.75\bar{\epsilon}_{0}$, and this is fit to experimental data on the
rat flexor hallucis brevis muscle in [63].

Fibres in the aponeurosis and tendon cannot be activated. The passive along-
fibre properties for the aponeurosis and tendon are the same, and are given by
a cubic polynomial fit to experimental data on the human gastrocnemius
reported in [44], see also [49].

$\hat{\sigma}_{\text{pass,apo/ten}}(\lambda)=\left\\{\begin{aligned}
&0&&0\leq\lambda\leq 1.0\\\
&515.882034(\lambda-1.0)^{2}+0.01(\lambda-1.0)+0.01&&1.0\leq\lambda\leq
1.01\\\
&600.590242(\lambda-1.01)^{2}+0.327640(\lambda-1.01)+0.06168820&&1.01\leq\lambda\leq
1.02\\\
&-9.975321(\lambda-1.02)^{2}+22.3394455(\lambda-1.02)+0.2250236&&1.02\leq\lambda\leq
1.15\\\ &19.7458618(\lambda-1.15)+2.960568&&\lambda>1.15,\end{aligned}\right.$
(30)