A continuum model of skeletal muscle tissue with loss of activation
Giulia Giantesio Alessandro Musesti


Low Rank Magnetic Resonance Fingerprinting
    Yonina C. Eldar
    
==========================================




We present a continuum model for the mechanical behavior of
  the skeletal muscle tissue when its functionality is reduced due to
  aging. The loss of ability of activating is typical of the geriatric
  syndrome called sarcopenia.  The material is described by a
  hyperelastic, polyconvex, transverse isotropic strain energy
  function. The three material parameters appearing in the energy are
  fitted by least square optimization on experimental
  data, while incompressibility is assumed through a Lagrange multiplier
  representing the hydrostatic pressure. The activation of the muscle
  fibers, which is related to the contraction of the sarcomere, is
  modeled by the so called active strain approach.  The loss of
  performance of an elder muscle is then obtained by lowering of some
  percentage the active part of the stress. The model is implemented
  numerically and the obtained results are discussed and graphically
  represented.




§ INTRODUCTION


Skeletal muscle tissue is one of the main components of the human
body, being about 40% of its total mass. Its principal role is the
production of force, which supports the body and becomes
movement by acting on bones. The mechanism by which a muscle
produces force is called activation.

Skeletal muscle tissue is a highly ordered hierarchical structure. The
cells of the tissue are the muscular fibers, having a length up to
several centimeters; they are organized in fascicles, where every
fiber is multiply connected to nerve axons, which drive the activation
of the tissue. Connective tissue, which is essentially isotropic,
fills the spaces among the fibers. Every fiber contains a
concatenation of millions of sarcomeres, which are the fundamental unit
of the muscle. With a length of some micrometers, a sarcomere is
composed by chains of proteins, mainly actin and myosin, which can
slide on each other. This sliding movement produces the contraction of
the sarcomere and, ultimately, the contraction of the whole muscle and
the production of force and movement.

The aim of this Chapter is to propose a mathematical model of skeletal
muscle tissue with a reduced activation, which is typical of a
geriatric syndrome named sarcopeniasarcopenia <cit.>.  About
thirty years ago, the term sarcopenia (from Greek sarx or flesh and
penia


or loss) has been introduced in order to describe the age-related
decrease of muscle mass and performance. Sarcopenia has since then
been defined as the loss of skeletal muscle mass and strength that
occurs with advancing age, which in turn affects balance, gait and
overall ability to perform even the simple tasks of daily living such
as rising from a chair or climbing steps. According to
<cit.>, sarcopenia affects more than 50 millions
people today and it will affect more than 200 millions in the next 40
years. There is still no generally accepted test for its diagnosis and
many efforts are made nowadays by the medical community to better
understand this syndrome. Therefore it is desirable to build a
mathematical model of muscle tissue affected by sarcopenia.  However,
to the best of our knowledge, in the biomathematical literature the
topic of loss of activation has never been addressed.

In order to use the valuable tools of Continuum Mechanics, during the
last decades the skeletal muscle tissue has been often modeled as a
continuum material <cit.>, 
which is usually assumed to be transversely
isotropic and incompressible. The former assumption is motivated by
the alignment of the muscular fibers, while the latter is ensured by
the high water content of the tissue (about 75% of the total volume).
Moreover, in view of some experimental tests, the material is assumed
to be nonlinear and viscoelastic. Focusing our attention only on the
steady properties of the tissue, here we neglect the viscous effects
and we set in the framework of hyperelasticity.

In the model that we propose, there are three constitutive
prescriptions: one for the hyperelastic energy when the tissue is not
active (passive energy), one for the activation and one for the loss
of performance. As far as the passive part is concerned, we assume an
exponential stress response of the material, which is customary in
biological tissues. The particular form that we choose, being a slight
simplification of the one proposed in <cit.>, has the 
advantage of being polyconvex and coercive, giving mathematical
soundness to the model and stability to the numerical
simulations.

A recent and very promising way to describe the activation is the
active strain approach, where the extra energy produced by the
activation mechanism is encoded in a multiplicative decomposition of
the deformation gradient in an elastic and an active part (see
Section <ref>). Unlike the classical
active stress approach, in which the active part of the stress
is modeled in a pure phenomenological way and a new term has to be
added to the passive energy, the active strain method does not change
the form of the elastic energy, keeping in particular all its
mathematical properties. Moreover, at least in the case of skeletal
muscles, the active strain approach seems to be more adherent to the
physiology of the tissue, in the sense that at the molecular level the
production of force is actually given by a deformation of the
material, thanks to the contraction of the sarcomeres. The active part
of the deformation gradient is a mathematical representation of such a
contraction. The multiplicative decomposition of the deformation
gradient has been applied to an active striated muscle in
<cit.>. However, this decomposition involves only a
part of the whole elastic energy, which is written as the sum of two
terms for the case of a fiber-reinforced material. As far as we know,
the active strain approach has never been previously applied to the
whole elastic energy of a skeletal muscle tissue.  As a drawback, the
active strain approach can be a source of some technical difficulties;
for instance, in our case fitting the model on the experimental data
is not so simple, see eq. (<ref>).

Furthermore, we consider the loss of performance, which is one of the
novelties of our model. Unfortunately, there are no experimental data
on the elastic properties of a sarcopenic muscle tissue, at least to
our knowledge; hence we adopt the naive strategy of reducing the
active part of the stress (which is the difference between the stress
of the material with and without activation) by a given percentage, represented
by the damage parameter d (see
Section <ref>). In this way, there is a single
parameter in the model which concisely accounts for any effect of the
disease.

The proposed model can be numerically implemented using finite element
methods. In Section <ref> we present some results
obtained using FEniCS, an open source collection of Python
libraries. Actually, we consider a cylindrical geometry with radial
symmetry, so that the numerical domain is two-dimensional and the
computational cost is reduced. As far as the boundary conditions are
concerned, we prescribe the displacement on the bases of the cylinder
and let the lateral surface traction-free. Such simulations show that
the experimental results of <cit.> on the passive and active
stress-strain healthy curves, obtained in vivo from a tetanized
tibialis anterior of a rat, can be well reproduced by our
model. Further, the behavior of the tissue when d increases is
analyzed. An ongoing task is to perform a finite element
implementation of the model when generic loads are applied, and to
consider a realistic three-dimensional muscle mesh.  We are now
developing a truly hyperelastic model, where the expression of the
stress takes into account also the dependence of the activation on the
deformation gradient. Actually, in this chapter the stress is computed
as the derivative of the hyperelastic energy keeping the active part
of the deformation gradient fixed.

In the future, it will be very interesting to find some connections
between the damage percentage (the parameter d) and other
physiological quantities, such as the mass of the muscular tissue or
the neuronal activity. Another important topic will be the application
of some homogenization techniques in order to deduce an improved
constitutive equation for the skeletal muscle starting from its
microstructure.






§ CONSTITUTIVE MODEL


Skeletal muscle tissue is characterized by densely packed muscle
fibers,muscle fibers which are arranged in fascicles. Filling the spaces between
the fibers and fascicles, connective tissue surrounds the muscle and
it is responsible of the elastic recoil of the muscle to elongation.
Besides a large amount of water, the fibers themselves contain titin,
actin and myosin filaments. The latter two sliding elements form the
actual contractile component of the muscle, which is called
sarcomere. Since the fibers locally follow a predominant
unidirectional alignment, transverse isotropy with respect to that
main direction can be assumed. We hence begin by modelling the
skeletal muscle tissue as a transversely isotropic nonlinear
hyperelastic material with principal direction , which follows
the alignment of the muscle fibers.



 §.§ Passive model

passive energy
Let  denote the deformation gradient tensor, =^T the right
Cauchy-Green tensor and =⊗ the so called structural
  tensor.structural tensor  If Ω denotes the reference configuration occupied
by the muscle, we describe its passive behavior by choosing a
hyperelastic strain energy function

    ∫_Ω W()dV,

where the strain energy density is of the form 

    W()=μ/4{1/α[e^α(I_p-1)-1]+K_p-1},

with

    I_p=w_0/3()+(1-w_0)(),  K_p=w_0/3(^-1)+(1-w_0)(^-1).

Here μ is an elastic parameter and α and w_0 are positive
dimensionless material parameters. The generalized invariants I_p
and K_p are given by a weighted combination of the isotropic and
anisotropic components; in particular, w_0 measures the ratio of
isotropic tissue constituents and 1-w_0 that of muscle
fibers. Moreover, the term () represents the squared
stretch in the direction of the muscle fiber and is thus associated with
longitudinal fiber properties, while the term (^-1)
describes the change of the squared cross-sectional area of a surface
element which is normal to the direction  in the
reference configuration and thus relates to the transverse
behavior of the material <cit.> (see Fig. <ref>).  

One of the mathematical features of the energy density
(<ref>) is that it is polyconvex and coercive
<cit.>, hence the equilibrium problem with mixed boundary
conditions is well posed.

We remark that  is the identity tensor  in the reference
configuration, so that I_p=K_p=1, i.e. we have the energy- and
stress-free state of the passive muscle tissue (see <cit.>).

The high content of water is responsible of the nearly incompressible
behavior which is experimentally reported for muscle fibers, so that
we can assumeincompressibility

    = 1.


As is customary in hyperelasticity, the first
Piola-KirchhoffPiola-Kirchhoff stress tensor stress
tensor, known as nominal stress tensor, can be directly
computed by differentiating the strain energy function:

    =   W-p^-T=2W-p^-T= 
       =   μ/2{e^α(I_p-1)[w_0/3+(1-w_0)]-^-1[w_0/3+(1-w_0)]^-1}-p^-T,

where p is a Lagrange multiplier associated with the hydrostatic
pressure which results from the incompressibility constraint
(<ref>).

The material parameters of the model can be obtained from real data.
More precisely, concerning the elastic parameter μ, we use the
value given in <cit.>, while the other two parameters have been
obtained by least squares optimization using the experimental data by
Hawkins and Bey <cit.> about the stretch response of a
tetanized tibialis anterior of a rat (see
Fig. <ref>). In Table <ref> we furnish the
values of the parameters.



 

 


We remark that the strain energy function (<ref>) is a
slight simplification of the one proposed by Ehret, Böl and Itskov in <cit.>:

    W_EBI()=μ/4{1/α[e^α(I_p-1)-1]+
    1/β[e^β(K_p-1)-1]},

where α=19.69, β=1.190, w_0=0.7388.  Actually, our
simplification consists in linearizing the term related to K_p,
which describes the transverse behavior. This is motivated by the
fact that the parameter β is much smaller than α. In
Fig. <ref> we can see the comparison between the nominal
stress in the direction of the stretch of the two models when the
muscle fibers are elongated in their direction.
 



 §.§ Active model

activationOne of the main features of the skeletal muscle tissue is its ability
of being voluntarily activated. Skeletal muscles are activated
through electrical impulses from motor nerves; the activation triggers
a chemical reaction between the actin and myosin filaments which
produces a sliding of the molecular chains, causing a contraction of
the muscle fibers.

During the last decades, many authors tried to mathematically model
the process of activation, mainly with two different approaches (for a
review see <cit.>).  The most famous approach followed in the
literature is called active stress and it consists in adding an
extra term to the stress, which accounts for the contribution
given by the activation (see for example
<cit.>).  However, this is an ad
  hoc method, usually not related to the sliding movement of the
filaments in the sarcomeres, which is the main mechanism of
contraction at the mesoscale.

More recently, the active strainactivation!active strain 
approach was proposed by Taber
and Perucchio <cit.> in order to describe the activation of
the cardiac tissue, following previous theories of growth and
morphogenesis, as well as several models of plasticity. The method for
soft living tissues is explained in <cit.>.
Differently from the active stress approach, this method does not
change the form of the strain energy function; rather, it assumes that
only a part of the deformation gradient, obtained by a multiplicative
decomposition, is responsible for the store of elastic energy.  This
method is related to the biological meaning of activation and can be
reasonably adopted also in our case. To the best of our knowledge, the
active strain approach has never been followed for the skeletal muscle
tissue in literature.

We begin by rewriting the deformation gradient as =_e_a, where
_e is the elastic part and _a describes the active
contribution (see Fig. <ref>). 
 
The active strain _a represents a change of the reference
volume elements due to the contraction of the sarcomeres, so that it
does not contribute to the elastic energy.  


A reference volume element, distorted by _a, needs a further
deformation _e to match the actual volume element,
which accommodates
both the external forces and the active contraction. Notice that neither _a nor
_e need to be the gradients of some displacement, that is, it is
not necessary that they fulfill the compatibility condition
curl_a=0 or curl_e=0.

The volume elements are modified by the internal active forces without
changing the elastic energy, hence the strain energy function of the
activated material has to be computed using _e=_e^T_e and
taking into account _e=_a^-1. If _a=χ_a
for some displacement χ_a, then from
Fig. <ref> by a change of variables it is easy to see that

    ∫_χ_a(Ω)W(_e)dV=∫_ΩW(^-T_a_a^-1)(_a)dV.

The right-hand side of the previous equation is well defined also when
_a does not come from a global displacement, and it describes
the strain energy of the active body.
We then obtain the modified hyperelastic energy density  

    W() = (_a)W(_e) = 
    (_a)W(^-T_a_a^-1).


We now have to model the active part _a. 
Since the activation of the muscle consists in a contraction along the fibers, we choose

    _a=-γ⊗,

where 0≤γ <1  is a dimensionless parameter representing
the relative contraction of activated fibers (γ=0 meaning no activation).
Then the modified strain energy density becomes 

    W()=(1-γ)W(_e)=(1-γ)μ/4{1/α[e^α(I_e-1)-1]+K_e-1
    }, 
          I_e=w_0/3(_e)+(1-w_0)(_e),   
    K_e=w_0/3(_e^-1)+(1-w_0)(_e^-1).

The corresponding first Piola-Kirchhoff stress tensorPiola-Kirchhoff stress tensor is given by

    =      _aW_e_a^-1-p^-T= 
    
     =      μ/2(1-γ)_e{e^α(I_e-1)[w_0/3+(1-w_0)]-_e^-1[w_0/3+(1-w_0)]
     _e^-1}_a^-1
          -p^-T,

where p accounts for the incompressibility constraint
=1. Notice that, since the activation (<ref>) does not
preserve volume and the material has to be globally incompressible,
one has that _e≠ 1, so that the material is elastically
compressible. As far as the strain energy density is concerned, a
factor (1-γ) appears in (<ref>) which keeps into
account the compressibility of _a. It would be interesting to
study also other kinds of passive energies, involving the quantity
, in order to better describe the elastic compressibility of
the material.

In Fig. <ref> we represent, for several values of the parameter
γ, the stress-strain curve for a uniaxial tension
along the fibers. If the muscle is activated (γ >0), then (the
absolute value of) the stress increases with γ and the value of
the stretch such that the stress is zero becomes less than one.

 



§ MODELLING THE ACTIVATION ON EXPERIMENTAL DATA

activation

The activation parameter γ, which was assumed constant in the
previous section, in fact usually depends on the deformation
gradient. In typical experiments on a tetanized skeletal muscle it is
apparent that the contraction of the fibers due to activation varies
with their stretch, reaching a maximum value and then
decreasing. Fig. <ref> shows the qualitative
relation between the elongation and the developed stress.
This section will be devoted to taking into account this
phenomenon. Specifically, the expression of γ will be
determined matching an experiment-based relation between stress and
strain with our model (<ref>).

In order to find the relation between stress and strain, the
experiments in vivo are usually performed in two steps. First,
the stress-strain curve is obtained without any activation
(passive curve). Second, by an electrical stimulus the muscle
is isometrically kept in a tetanized state and the total
  stress-strain curve is plotted. The last curve, which is
qualitatively represented in Fig. <ref>, depends on
the reciprocal position of actin and myosin chains. By taking the
difference of the two curves one can obtain the active curve,
describing the amount of stress due to activation. It is useful to
find a mathematical expression of such a curve, in order to take into
account the experimental behavior of the active contraction. This
issue has already been addressed in several papers, see
e.g. <cit.>.

 

Denoting with λ the ratio between the current length of the
muscle and its original length, we assume the active curve to be of
the form

    P_act(λ) =
    {   P_optexp[-k(λ^2-λ_opt^2)^2/λ-λ_min]
       if λ>λ_min,
    
       0    otherwise,
    .

where λ_min is the minimum stretch value after which the
activation starts (i.e.  the lower bound for the stretch at which the
myofilaments begin to overlap) and k is merely a fitting
parameter. The coordinates (λ_opt,P_opt) identify the
position of the maximum of the curve.  As it is explained in
<cit.>, the value of P_opt takes into account some
information at the mesoscale level, such as the number of activated
motor units and the interstimulus interval; according to the
literature <cit.>, it is set at P_opt=73.52
kPa. The numerical values of the other three parameters, deduced
through least squares optimization on the data reported in
<cit.>, are given in Table <ref>.

The expression (<ref>) has the advantage of describing the
asymmetry between the ascending and descending branches of the active
curve obtained in <cit.>. Indeed, even if the asymmetry is not
so evident in their curve, due to the fact that there are only few
data on the descending branch, it is a typical feature of several
experimentally measured sarcomere length-force relation.  Moreover, as
one can easily see in Fig. <ref>, the convex behavior of the
data nearby λ_min is well fitted.
 




 §.§ The activation parameter γ as a function of the elongation

activation!parameter of

Now our aim is to obtain P_act(λ) given in (<ref>)
from the model described in Section
<ref>. In order to reach our purpose, we
have to model the activation parameter γ as a function of the
stretch.

As in the experiments of Hawkins and Bey <cit.>, let us consider a uniaxial simple tension along the fibers. For simplicity, we assume that the fibers follow the direction 
=_1. 
Since the skeletal muscle tissue is modeled as an incompressible transversely isotropic material, the general form of the 
deformation gradient  is given by  

    =
    [      λ      0      0;      0 1/√(λ)      0;      0      0 1/√(λ) ].


Then using the notation introduced in Section <ref>, one has

    _e   =
    [ λ^2/(1-γ)^2           0           0;           0         1/λ           0;           0           0         1/λ ], 
       I_e   =w_0/3[λ^2/(1-γ)^2+2/λ]+(1-w_0)λ^2/(1-γ)^2, 
       K_e   =w_0/3[(1-γ)^2/λ^2+2λ]+(1-w_0)(1-γ)^2/λ^2.


In this case, it is convenient to look at the strain energy as a function of the stretch λ and the activation parameter γ:

    W(λ,γ)=(1-γ)W(λ,γ)=(1-γ)μ/4{1/α[e^α(I_e-1)-1]+K_e-1
    }.

 
Then the nominal stress along the fiber direction is given by

    P_tot(λ,γ):= Wλ
     =(1-γ)μ/4[I_e'
    e^α(I_e-1) +K_e'],

where 

    I_e'   =I_eλ=2w_0/3[λ/(1-γ)^2-1/λ^2]+2(1-w_0)λ/(1-γ)^2,
       K_e'   =K_eλ=2w_0/3[-(1-γ)^2/λ^3+1]
    -2(1-w_0)(1-γ)^2/λ^3.


We can get the passive stress by setting γ=0: 

    P_pas(λ):=P_tot(λ,0)
    =   μ2   {[(1-2/3w_0)λ-w_0/31/λ^2]e^α[(1-2/3w_0)λ^2+w_0/32/λ-1].
       -   .(1-2/3w_0)1/λ^3+w_0/3}.

We remark that the values of P_tot and P_pas can also be
obtained by computing the first component of the stress given by
(<ref>) and (<ref>) after finding the
hydrostatic pressure from the conditions
P_22=P_33=P_22=P_33=0 (traction-free
lateral surface).

Our aim is to find the value of γ such that

    P_tot(λ, γ)= P_act(λ)+P_pas(λ),

where P_act (λ) is given by (<ref>). Unfortunately,
this leads to an equation for γ which cannot be explicitly
solved:

    (1-γ){[(1-2/3w_0)λ/(1-γ)^2
    -w_0/31/λ^2]
        e^α[
        (1-2/3w_0)λ^2/(1-γ)^2+w_0/32/λ
        -1].
    .+w_0/3
    -(1-2/3w_0)(1-γ)^2/λ^3}
        =2/μ[P_act(λ)+P_pas(λ)].

However one can employ standard numerical methods and plot the
solution. Fig. <ref>_1, which is obtain by a bisection
method, shows γ as a function of λ.
We remark that γ vanishes before λ_min, indeed in
this region there is no difference between total and passive
stress. The corresponding behavior of the stresses is plotted in
Fig. <ref>_2, which is very similar to the
representative plot of Fig. <ref>.

We emphasize that the previous model is not strictly
hyperelastic, since in the expression of the
stress (<ref>) the derivative of γ with
respect to 𝐅 has been neglected. We are now working on a
truly hyperelastic model, which can be useful for some numerical
implementations.

 



 §.§ Loss of activation
 
activation!loss of

We now want to describe from a mathematical point of view the loss of
performance of a skeletal muscle tissue. As we have already explained
in the Introduction, this is one of the main effects of
sarcopenia,sarcopenia which is a typical syndrome of
advanced age.

In <cit.> it is remarked that aging is
associated with changes in muscle mass, composition, activation and
material properties. In sarcopenic muscle, there is a loss of motor
units via denervation and a net conversion in slow fibers, with a
resulting loss in muscle power. Hence, the loss of performance of a
sarcopenic muscle can be described as a weakening of the activation of
the fibers.

Unfortunately, as far as we know, there are no experimental data
describing a uniaxial simple tension along the fibers of a sarcopenic
muscle. For this reason, we try to describe the loss of activation by
a parameter d which lowers the curve P_act(λ) given by
(<ref>). The parameter d describes the percentage of
disease or damage: if d=0, then the muscle is
healthy. In order to get our aim, we multiply the function
P_act(λ) by the factor 1-d, as one can see in
Fig. <ref>. Notice that such a choice can be overly simple:
for instance, it implies that the maximum is always attained at
λ_opt, even if there is no experimental evidence of
that. However, the presence of d allows to describe, at least
qualitatively, the loss of performance of a muscle, which is one of the
goals of our model.






§ NUMERICAL VALIDATION



Finally, we simulate numerically the contraction and the elongation of
a slab of skeletal muscle tissue represented by a cylinder. We assume
radial symmetry, so that the mesh is a rectangle.  The ends of the
cylinder are assumed to remain perpendicular to the axial direction.
The rectangle is modeled by the hyperelastic model presented in the
previous sections. The active contractile fibers are aligned along the
length of the rectangle, which coincides with _1.  The passive
and active material parameters are given in Tables <ref> and
<ref>, respectively. Concerning the boundary conditions, the
cylinder is fixed at one end and elongated to a given length, in order
to recreate the situation of the experiments reported in
<cit.>. The lateral surface is assumed to be tension-free.

The analysis is performed by using the computing environment
FEniCS. The FEniCS ProjectFEniCS <cit.> is a collection of
numerical software, supported by a set of novel algorithms and
techniques, aimed at the automated solution of differential equations
using finite element methods.

As it is explained in Section <ref>, one of
the main features of our model is the dependence of the activation
parameter γ on the stretch λ. The function
γ(λ) solves the implicit equation
(<ref>), which ensures that the corresponding
stress curves fit the experimental data. However, even if this
equation can be solved using numerical methods, it is interesting to
find an explicit function in order to analyze qualitatively the active
model and to run the simulations in FEniCS. Moreover, the explicit
function γ(λ) has to be very precise, since a slight
error on γ deeply affects the behavior of the total stress.
Hence, it is reasonable to relate the expression of γ to 
the material parameters and the quantities involved in
(<ref>).
An idea is to isolate the exponential in (<ref>) 
and to express its exponent by a first step approximation of a
fixed-point method. We then obtain the following expression of
γ:

    γ(λ)=   {   a[√(1-2/3w_0/g(λ_min)/α+w_0/3λ_min)λ_min
    -√(1-2/3w_0/g(λ)/α+w_0/3λ)λ]
       if λ>λ_min,
    
       0    otherwise,
    .
     
       g(λ)=   lnα+α(1-w_0/λ)-1/2ln(1-2/3w_0/1/α+w_0/3λ)
          +ln{b2/μ[P_act(λ)+P_pas(λ)]+
    √(1-2/3w_0/1/α+w_0/3λ)[
    (1-2/3w_0)^2/1/α+w_0/3λ-λw_0/3]},

where a and b are dimensionless fitting parameters: a is related to the magnitude
of γ, while b acts on the curves (<ref>) and
(<ref>), which are the terms of the equation not depending on
γ.  Performing a least square optimization on the resulting
P_act, one gets a=1.0133 and b=0.2050.

Fig. <ref> shows the plot of the function
γ(λ) given in (<ref>) in comparison to the
numerical solution of equation (<ref>) obtained by
a bisection method.
 
Notice that the function defined in (<ref>) is continuous; in
particular we impose γ(λ_min)=0, so that the starting
value of activation does not change. Moreover, the function
approximates very well the numerical values of γ in the range
0.7<λ<1.5. However, the fitting is not so good when λ
becomes larger: for instance, the function is negative for
λ≥ 1.6. Nevertheless, the latter behavior of γ does not
influence too much the curve P_tot, since in that region
P_pas≫ P_act. Indeed, one can even neglect the activation for
large stretches.
The total stress response  is plotted in Fig. <ref>
in comparison to the data given in <cit.>.
 

Finally, it is interesting to run the simulations in the case of loss
of activation, i.e. when the damage parameter d varies. In
order to find the suitable activation function γ(λ), it
is sufficient to multiply the term P_act in (<ref>) by
(1-d).  As one would expect from Fig. <ref>, we have that
when d increases the activation γ decreases
(Fig. <ref>_1). This means that lowering the curve of
P_act results in a decrease of γ(λ), which leads to a
lowered total stress response. As one can see in
Fig. <ref>_2, the damage parameter mainly affects the
value of the stress in the region near λ_opt, where the
active stress reaches its maximum. However, the qualitative behavior of the stress curve
does not change, at least for d≤ 0.5. In particular, after a
plateau, the stress follows the exponential growth of the passive curve.
 



§ ACKNOWLEDGEMENT

  This work has been supported by the project Active Ageing and
    Healthy Living <cit.> of the Università Cattolica
  del Sacro Cuore and partially supported by GNFM (Gruppo
  Nazionale per la Fisica Matematica) of INdAM (Istituto Nazionale di
  Alta Matematica).

The authors wish to thank the anonymous referees for their useful comments.


10

D:fenics
M. S. Alnæs, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg,
  C. Richardson, J. Ring, M. E. Rognes, and G. N. Wells.
The FEniCS Project Version 1.5.
Archive of Numerical Software, 100:9–23, 2015.

D:asasb
D. Ambrosi and S. Pezzuto.
Active Stress vs. Active Strain in Mechanobiology:
  Constitutive Issues.
Journal of Elasticity, 107:199–212, 2012.

D:blemker
S. S. Blemker, P. M. Pinsky, and S. L. Delp.
A 3D model of muscle reveals the causes of nonuniform strains in
  the biceps brachii.
Journal of Biomechanics, 38:657–665, 2005.

D:bol
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:489–504, 2008.

D:review
G. Chagnon, M. Rebouah, and D. Favier.
Hyperelastic Energy Densities for Soft Biological
  Tissues: A Review.
Journal of Elasticity, 120:129–160, 2015.

D:reportSarcopenia
A. J. Cruz-Jentoft, J. P. Baeyens, J. M. Bauer, Y. Boirie, T. Cederholm,
  F. Landi, F. C. Martin, J. P. Michel, Y. Rolland, S. M. Schneider,
  E. Topinková, M. Vandewoude, and M. Zamboni.
Sarcopenia: European consensus on definition and diagnosis.
Age and Ageing, 39:412–423, 2010.

D:ebi
A. E. Ehret, M. Böl, and M. Itskov.
A continuum constitutive model for the active behaviour of skeletal
  muscle.
Journal of the Mechanics and Physics of Solids, 59:625–636,
  2011.

D:ebipcb
A. E. Ehret and M. Itskov.
A polyconvex hyperelastic model for fiber-reinforced materials in
  application to soft tissues.
Journal of Materials Science, 42:8853–8863, 2007.

D:ei2009
A. E. Ehret and M. Itskov.
Modeling of anisotropic softening phenomena: Application to soft
  biological tissues.
International Journal of Plasticity, 25:901–919, 2009.

D:datib
D. Hawkins and M. Bey.
A Comprehensive Approach for Studying Muscle-Tendon
  Mechanics.
ASME Journal of Biomechanical Engineering, 116:51–55, 1994.

D:thomas
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:1–14, 2014.

D:hernandez
B. Hernández-Gascón, J. Grasa, B. Calvo, and J. F. Rodríguez.
A 3D electro-mechanical continuum model for simulating skeletal
  muscle contraction.
Journal of Theoretical Biology, 335:108–118, 2013.

D:Johansson
T. Johansson, P. Meier, and R. Blickhan.
A Finite-Element Model for the Mechanical Analysis of
  Skeletal Muscles.
Journal of Theoretical Biology, 206:131–149, 2000.

D:Lang2010
T. Lang, T. Streeper, P. Cawthon, K. Baldwin, D. R. Taaffe, and T. B. Harris.
Sarcopenia: etiology, clinical consequences, intervention, and
  assessment.
Osteoporos Int, 21:543–559, 2010.

D:martins
J. A. C. Martins, E. B. Pires, R. Salvado, and P. B. Dinis.
A numerical model of passive and active behavior of skeletal muscles.
Computer Methods in Applied Mechanics and Engineering,
  151:419–433, 1998.

D:giulio
A. Musesti, G. G. Giusteri, and A. Marzocchi.
Predicting Ageing: On the Mathematical Modelization of
  Ageing Muscle Tissue.
In G. Riva et al., editor, Active Ageing and Healthy Living.
  IOS press, 2014.
Chapter 17.

D:nardinocchiteresi
P. Nardinocchi and L. Teresi.
On the Active Response of Soft Living Tissues.
Journal of Elasticity, 88:27–39, 2007.

D:PTD
C. Paetsch, B. A. Trimmer, and A. Dorfmann.
A constitutive model for active-passive transition of muscle fibers.
International Journal of Non-Linear Mechanics, 47:377–387,
  2012.

D:progettob
G. Riva, P. Ajmone Marsan, and C. Grassi.
Active Ageing and Healthy Living.
IOS press, 2014.

D:sch
J. Schröder and P. Neff.
Invariant formulation of hyperelastic transverse isotropy based on
  polyconvex free energy functions.
International Journal of Solids and Structures, 40:401–445,
  2003.

D:taber
L. A. Taber and R. Perucchio.
Modeling Heart Development.
Journal of Elasticity, 61:165–197, 2000.

D:vanLeeuwen1991
J. L. van Leeuwen.
Optimum power output and structural design of sarcomeres.
Journal of Theoretical Biology, 149:229–256, 1991.

D:vanLeeuwen1992
J. L. van Leeuwen.
Muscle function in locomotion.
In Advances in Comparative and Environmental Physiology.
  Springer Heidelberg Berlin, 1992.
Chapter 7.

D:overviewSarco
S. von Haehling, J. E. Morley, and S. D. Anker.
An overview of sarcopenia: facts and numbers on prevalence and
  clinical impact.
Journal of Cachexia, Sarcopenia and Muscle, 1:129–133, 2010.



