An Analysis of Including the Evolution Law for the Serial Element in the Musculoskeletal Modelling

(1)

Biomedical Engineering

Master’s Thesis

An Analysis of Including the Evolution Law for the

Serial Element in the Musculoskeletal Modelling

Alexandra Roser

Supervisor:

Jonas St˚

alhand, PhD

IEI, Link¨oping University

Examiner:

Joakim Holmberg, PhD

IEI, Link¨oping University

LIU-IEI-TEK-A--19/03583--SE

Division of Solid Mechanics

Department of Management and Engineering (IEI)

The Institute of Technology

Link¨oping University

SE-581 83 Link¨oping, Sweden

(2)

Acknowledgements

Of all the people who have supported me throughout this thesis and my studies as a whole, I would like to show my appreciation to a few in particular.

First and foremost, I would like to thank my supervisor Jonas St˚alhand and examiner Joakim Holmberg for their endless effort and support. Your guidance and feedback is very much appreciated and I greatly enjoyed my thesis topic and my time at the division. I would also like to thank you for the countless interesting conversations and discussions, on as well as off topic.

A special thanks goes to my office neighbours Jan-Lucas and Christian for the excellent coffee and the effective training that brought me a lot of joy (and pain). I would like to thank Jan-Lucas in particular for being my debugging rubber duck in times of need.

Last, but certainly not least, I would like to thank my family for the years of unconditional support.

Link¨oping, December 2019 Alexandra Roser

(3)

Nomenclature

Symbol Description Unit

C1 Specific muscle material parameter −

C2 Specific muscle material parameter −

F Externally applied force N

FD Damping force N

FM Muscle force N

FM

0 Maximum isometric force N

FM

max Maximum muscle force N

FT _{Tendon force} _N

FP E _{Parallel element force} _N

Fm Muscle force in the AnyBody Modeling System N

Ft Tendon force in the AnyBody Modeling System N

Fpe Parallel element force in the AnyBody Modeling System N

JT Tendon shape parameter −

LM Muscle length m

LM

B Muscle belly length m

LM

max Maximum muscle length m

LM

min Minimum muscle length m

LM

opt Optimal muscle length m

LT _{Tendon length} _m

LT₀ Tendon length at maximum isometric force m

LTs Tendon slack length m

LT_isot Tendon length during isotonic contraction m

Lm Muscle length in the AnyBody Modeling System m

Lt Tendon length in the AnyBody Modeling System m

Lmt Total length in the AnyBody Modeling System m

L Total muscle length m

Mt Torque equalling ATP consumption to generate muscle force N m

P Power of the modified model N m/s

(4)

Strength Fully activated muscle force in the AnyBody Modeling System N

VM _{Muscle volume} _ml

W Weight N

Ψ Strain-energy N m

ΨT Strain-energy in the serial spring N m

ΨP E Strain-energy in the parallel spring N m

α Muscle specific constant N

β Muscle specific constant m/s

˙ LM

max Maximum contraction speed m/s

γ Pennation angle of muscle fibres ◦

ˆ

R2 _{Modified coefficient of determination} ₋

ˆ

κ Muscle specific parameter −

κ Constant −

λM Muscle stretch −

µ Mean of Gaussian force-length dependence −

ω Angular velocity of the friction clutch 1/s

φ Force dependence on activated cross-bridges −

σ Standard deviation of Gaussian force-length dependence −

τ Afterload N

τ0 Maximum contractile force at current muscle length N

εT

0 Tendon strain at maximum isometric force −

ζ Variable N s/m

a Activation variable −

fL(LM) Force-length dependence −

fv( ˙LM) Force-speed dependence −

k Constant −

r Radius of the friction clutch m

t Time s

Abbreviations Ca2+ _{Calcium ion}

ADP Adenosine Diphosphate ATP Adenosine Triphosphate FT Fast twitch muscle fibre type

(5)

PCSA Physiological cross-sectional area of the muscle ST Slow twitch muscle fibre type

Subscripts

m Muscle in the AnyBody Modeling System

mt Muscle and tendon in the AnyBody Modeling System pe Parallel element in the AnyBody Modeling System t Tendon in the AnyBody Modeling System

Superscripts CE Contractile element D Damper M Muscle P E Parallel element SE Serial element T Tendon iv

(6)

1 Introduction

The musculoskeletal system is composed of bones, skeletal muscles, joints, and connective tissue such as tendons and ligaments. Their interactions provide the body with stability and allow for movement, but also protect vital organs. Muscles are biological actuators that generate movement by contracting and relaxing. They are connected to the bones through tendons at both ends of the muscle. The location of the attachment to the bone determines the direction of the movement. While one single monoarticular muscle is only able to perform one motion, the interaction of muscle groups allows for complex movements.

In order to simulate muscle contraction, several mechanical models have been developed. One of the first models was proposed byHill(1938), and it is still used today. Hill performed a set of experiments on frog muscles and was able to relate tiny temperature changes in the muscle to the generated muscle forces. Based on a thermodynamical argument, he proposed a mechanical muscle model and an equation relating muscle force to its contraction speed. The muscle model consists of a contractile element (CE) and two springs: one in parallel and one in series, see Fig.1.1(left). The contractile element is usually taken to be the muscle itself in which the force is generated, while the serial element (SE) can be interpreted as the tendon, and the parallel element (PE) as passive resistance in the muscle itself, e.g. due to connective tissue. Strictly speaking,Hill(1938) only concluded that there is an elastic element in series with the contractile element but it is common practice to include a parallel element as well (e.g.Christensen et al.,2002;Siebert et al.,

2008).

Since the serial element can be arranged in two ways, the Hill model is not unique. A common model is given by arranging the serial element together with the contractile element in parallel to the parallel element (Winters, 1990), see Fig. 1.1 (right). It has, however, been shown (Siebert et al., 2008) that the ‘classic’ Hill model to the left represents the behaviour of skeletal muscles more accurately than the Hill model to the right, which is why the former will be the focus of this thesis. From now on, it will simply be referred to as the Hill model.

At this point, it should be noted that several different terms for the contractile, serial, and parallel elements exist and are used simultaneously throughout the literature. While the term ‘muscle’ commonly refers to the whole muscle including the tendons, these structures are kept separate herein. The contractile element, or the contractile unit, is the pure muscle part in which force is generated, and will be referred to as the muscle. The serial element is modelled as a spring and represents the tendon. The parallel element is also modelled as a spring and is sometimes referred to as passive element. In the following, parameters related to the contractile element will be denoted by a superscribed M to indicate muscle, while parameters related to the serial element are denoted by a superscribed T for tendon, and parameters for the parallel element are denoted

CE SE PE Force Force Muscle length Force Force CE SE PE Tendon length Muscle length Total length

Figure 1.1: Two different arrangements of the contractile element (CE) and the serial (SE) and parallel elements (PE) in the Hill model.

(9)

CHAPTER 1. INTRODUCTION 2

by a superscribed P E.

For the sake of completeness, it should be noted that in a physiological configuration, there is usually one tendon attached to either side of the muscle. For simplicity, these two tendons are commonly replaced by a single spring, as this replacement is mechanically equivalent.

As can be seen in the left Hill model in Fig. 1.1, the contractile and parallel elements are of equal length, which is referred to as the muscle length. Together with the tendon length it forms the total length. The speed at which a muscle contracts is given by the rate of change of its total length with respect to time.

Why is it important to know the muscle length and its change in time? In his experiments under isometric conditions,Hill(1938) found that the muscle force is dependent on both the muscle length and the speed at which it contracts. Only at the so-called optimal muscle length, the muscle force would reach its maximum isometric force. At lengths falling below a minimum length or exceeding a maximum length, the muscle is unable to generate any force at all.

The generated force also decreases with decreasing muscle speed and the muscle is unable to generate force at speeds exceeding the maximum contraction speed, where the contraction speed is defined negative. Additionally, the greatest possible muscle force during contraction is the maximum isometric force achieved at zero speed. Professional athletes, for example cyclists, make use of this knowledge when adjusting the seat height of their bicycles (influencing the length of the thigh muscles) and when choosing a gear for an optimal pedalling rate (speed dependence) to maximise the generated force (Herzog,2000).

In so-called inverse muscle dynamics, the inner forces in the muscle are computed from the external force and the body movement. The measurable quantities in muscle experiments are usually the total length, the total speed, and the tendon force. However, the split of the total length into the muscle length and tendon length is not unique. For an isometric case, during which the total length is constant, there is an infinite number of combinations of muscle and tendon length whose sum equals the constant total length. In order to achieve a unique solution for each given total length, the tendon length is assumed to be constant, which is why the tendon speed is mostly considered negligible in musculoskeletal simulation tools (Holmberg, 2012). This assumption is motivated by the high tendon stiffness which results in a possibly insignificant change in tendon length during contraction (Holmberg, 2012).

Although the tendon exhibits a high stiffness, its length does not stay constant during contrac-tion, which is why the aforementioned assumption introduces errors of unknown magnitude to the simulations. It is therefore important to assess whether the tendon length can safely be assumed to be constant during muscle contraction, and if so, how large the introduced errors are.

In this thesis, the standard Hill model is modified to a more accurate model in which the mus-cle speed is numerically calculated taking into account a non-zero tendon speed. The simulations are carried out in Matlab®(The Mathworks, Natick, MA, USA) for a simple biceps curl. The results are eventually compared to simulations done in the commercial musculoskeletal simulation software The AnyBody Modeling System (AnyBody Technology, Aalborg Ø, Denmark) in order to study the effect of assuming a negligible tendon speed.

The value of the Hill model lies in its ability to give an acceptable accuracy despite its simplicity. It should, however, be noted that several limitations are introduced by e.g. ignoring geometrical effects and assuming homogeneity and one-dimensionality. Furthermore, as additional shortcom-ings of the model, thermal effects are disregarded and muscle activation is only considered as an instantaneous on/off switch.

(10)

2 Muscle Physiology

To understand the form of the force-length dependence observed in muscles, the structure and the underlying biological and chemical processes in a muscle fibre have to be regarded. In this section, some basic muscle physiology is presented. The intention is not to give a complete picture but to introduce concepts which will help the reader to understand the modelling. For an in-depth treatment, the reader is referred to, e.g.,Tortora and Derrickson(2013).

2.1 Structure of a Muscle

A schematic picture of the anatomical structure of a skeletal muscle can be seen in Fig.2.1. It is composed of individual muscle fibres that run either in parallel to the force-generating axis or at an angle to it (pennation angle).

Each individual muscle fibre is a cell which is filled with cellular fluid, called the sarcoplasm, several nuclei, and cell organelles, e.g. the mitochondria. The mitochondria in the cell produce adenosine triphosphate (ATP), an energy-carrying molecule needed for contraction. The muscle fibres are enclosed by a cell membrane, the sarcolemma. The sarcoplasm in a muscle fibre is also filled with several myofibrils, which are surrounded by a netlike organelle, the sarcoplasmic reticulum. The sarcoplasmic reticulum stores calcium ions (Ca2+_{) and releases them into the}

sarcoplasm after the muscle cell has been stimulated by an action potential, provoking muscle contraction. Along their length, the myofibrils themselves are compartmentalised into so-called sarcomeres, which can be described as the fundamental functional unit of a muscle (McMahon,

1984, p. 55). The sarcomeres are separated by the Z-disc and consist of thick and thin filaments arranged in layers. The thick and thin filaments are protein structures that are composed of the contractile proteins myosin and actin, respectively (see Fig. 2.2). The myosin filaments are transversely connected to each other by a protein structure called the M-line. In the presence of Ca2+ _{and using ATP, the myosin heads of the thick filament can attach to binding sites on the}

actin molecules of the thin filament, forming so-called cross-bridges that enable muscle contraction.

Bone Tendon

Skeletal muscle Fascicle

Muscle fibre Sarcomere Sarcolemma Sarcoplasm Myofibril Nucleus Actin Myosin Z-disc Mitochondrion

Figure 2.1: Structure of a skeletal muscle (reproduced fromTortora and Derrickson,2013, pp. 330-332).

2.2 Sliding Filament Theory

In 1954,Huxley and Hanson (1954) and Huxley and Niedergerke(1954) discovered, unrelated to each other and independently of one another, that filament length did not change during muscle contraction, but the filaments were sliding past each other. As a result, they postulated the sliding

(11)

CHAPTER 2. MUSCLE PHYSIOLOGY 4

filament theory. Gordon et al.(1966) provided further evidence for the sliding filament theory and also proposed a hypothesis for the force-length dependence in skeletal muscles based on the sliding filament theory.

2.2.1 Muscle Contraction

An electrical nerve impulse is sent to the muscle over nerve pathways during voluntary contraction. At the end of those pathways, the impulse triggers the release of neurotransmitters, which in turn cause a transient reduction in the cell membrane potential (action potential) and an influx of Ca2+_{. As the muscle action potentials subsequently travel along the sarcolemma, calcium ions}

are released from the sarcoplasmic reticulum into the sarcoplasm. The intracellular calcium ions initiate a reaction that frees the binding sites on the actin, making it possible for the myosin heads of the thick filament to attach to the thin filaments and start the contraction cycle. During the contraction cycle, ATP binds to the myosin head and is hydrolysed into adenosine diphosphate (ADP) and a phosphate group. The energy that is released during this breakdown process is used to attach the myosin head to the binding site on the actin, forming a cross-bridge and releasing the phosphate group. Once the cross-bridge is formed, it rotates in direction of the centre of the sarcomere, a process during which the ADP is released and force is produced. During this so-called power stroke, the thin filaments start to slide as they are connected to the thick filaments by the cross-bridges. As a result of the sliding thin filaments, the Z-discs to which they are connected are being drawn closer to each other, effectively shortening the sarcomere and therefore the whole muscle. The cross-bridge is not dissolved unless a new ATP molecule binds to the myosin head. This contraction cycle continues as long as ATP is supplied and the Ca2+ _{concentration in the}

sarcoplasm is high enough.

2.2.2 Myosin Cycling

During the contraction cycle, only a fraction of the myosin heads are attached to the thin filaments at the same time, whereas the rest of the myosin heads remain deattached until ATP is hydrolysed. Since the myosin heads take turns in being attached, there will always be active cross-bridges pulling the thin filaments towards centre of sarcomere as long as the contraction cycle is in progress.

Due to the structure of the single sarcomeres, more and more myosin heads are able to connect as the thin filaments slide towards the centre, leading to a steadily increasing force. At its peak, the thin and thick filaments are fully overlapping, which generates a maximum amount of force as a maximum amount of cross-bridges are formed. However, as the contraction cycle continues, the thin filaments keep sliding past that optimal point until the thin filaments themselves start overlapping, leading to a decreasing force. Eventually, the thick filaments will collide with the Z-discs at either end of the sarcomere, causing the thick filaments to buckle and fewer cross-bridges being able to form (McMahon,1984, p. 67).

The amount of cross-bridges that are active at the same time relate directly to the amount of force produced. This applies both to the amount of cross-bridges that are able to form due to the position of the filaments to each other, and due to the present calcium concentration.

Z-disc M-line Z-disc

Actin

Myosin

Cross-bridge

Figure 2.2: Thick filaments (myosin) and thin filaments (actin) in a sarcomere (reproduced from

(12)

2.3 Experimental Observations

The following chapter will cover the most important observations that were made during early muscle experiments.

2.3.1 Force-Length Dependence

The muscle force is related to the muscle length and is described by the characteristic force-length dependence of each muscle. The force-length dependence is obtained in isometric experiments by bringing the muscle to a constant length and fully activating the muscle. The measured muscle force can then be plotted over the muscle length as illustrated in Fig.2.3. The force-length dependence is therefore a depiction of the maximum force a muscle can generate at a particular length.

FM LM F₀M LM_min LM opt LMmax FM Time F₂M F₁M LM 1 LM2 LM₁ LM 2

Figure 2.3: Creation of a characteristic force-length response of a muscle under isometric conditions (reproduced fromHerzog, 2000). The maximum isometric force FM

0 is generated at the optimal

muscle length LM

opt, while the muscle is unable to generate force at muscle lengths below LMminand

above LM max.

Based on the sliding filament theory described in section 2.2, the graph of the force-length dependence in Fig.2.3can be explained and divided into three main sections.

Gordon et al. (1966) assigned different filament stages to certain sections in the force-length dependence, as shown in Fig.2.4. The linear descent in force during section I is related to the increasing amount of myosin heads that are able to form cross-bridges as the actin filaments start sliding towards the centre of the sarcomere. During the so-called plateau region in section II, the filaments are overlapping in a way that allows all myosin heads to attach to the actin, enabling the muscle to generate maximum force. The steep ascent in force in section III is attributed to the buckling of the myosin filaments as they reach the Z-disc at either end of the sarcomere until no force generation is possible any longer.

F0M Sarcomere length [µm] 1.0 1.5 2.0 2.5 3.0 3.5 4.0 I II III I II III Myosin Actin FM Z-disc Z-disc

Figure 2.4: Force-length dependence based on the position of the filaments (reproduced from

(13)

The force-length dependence was motivated based on findings at sarcomere level but can be scaled to represent the whole muscle because of the self-similarity of the sarcomeres and since the muscle behaves like the sarcomeres it is composed of.

Figure 2.5 shows the force-length dependence of the muscle, tendon, and parallel element, respectively, plotted over the muscle length LM. The total force is defined to be the sum of the forces of the contractile element, FM, and the parallel element, FP E, and is kinematically assumed to be the tendon force FT. It should be noted that these are the available muscle forces at maximum activation of the muscle. During actual movements, the muscle force and tendon force are usually smaller but always lie within their respective available curve. The passive parallel element force is however always proportional to the amount of stretching in the muscle, and therefore independent of muscle activation.

F

LM

F₀M

LM

min LMopt LMmax

FP E FT

FM

Figure 2.5: Characteristic force-length dependence of available forces in the muscle, tendon, and parallel element at a pennation angle of zero.

2.3.2 Force-Speed Dependence

The force-speed dependence describes the relation of the muscle force to the speed at which it contracts. Figure2.6shows that the generated muscle force FM decreases with decreasing muscle speed ˙LM _{until the muscle is unable to produce force at speeds lower than the maximum contraction}

speed − ˙LM

max. Additionally, it can be seen that the greatest possible muscle force during contraction

is the maximum isometric force FM

0 at zero speed. − ˙LMmax FM max Contraction Elongation FM 0 ˙ LM FM 0

Figure 2.6: Speed dependence of muscle force (reproduced fromHerzog,2000). In contrast to the accepted convention in the literature, ˙LM is positive in extension.

2.3.3 Hill Equation

Although Hill did not know about the sliding filament theory, he proposed an accurate muscle model, as outlined in chapter1. Based on his experiments, he also connected the muscle force to

(14)

the muscle’s contraction speed through the so-called Hill equation (Hill,1938) ˙ LM = β F M 0 − FM FM _{+ α} , L˙ M _{< 0,} _(2.1)

where α and β are constants. In the original quick-release experiments, FM _{is the constant}

(15)

3 Muscle Models

3.1 Classical Skeletal Muscle Model

The standard model for muscle contraction in Fig.3.1 has been used in many studies, see e.g.,

G¨unther et al. (2007); Manal and Buchanan (2013); Kim et al. (2014). It has a high accuracy while remaining simple (Winters, 1990). Depending on the muscle in question, the contractile and parallel element can be angled at an angle γ to the serial element, representing the possible pennation of the individual muscle fibres. While some muscles have a greater pennation angle, others exhibit almost no pennation at all, resulting in a muscle fibre arrangement nearly parallel to the tendon. Due to their structure, more fibres can work in parallel at a higher pennation angle, allowing pennate muscles to generate greater muscle forces (McMahon,1984, p. 55). In this thesis, pennation of the muscle will be disregarded, i.e. γ = 0◦.

F T LT LM F M PE γ L

Figure 3.1: Classical Hill model with pennation.

3.1.1 Kinematics

As can be seen in the Hill model (Fig.3.1), the contractile and parallel elements are of equal length, which will be referred to as the muscle length LM_{. Together with the tendon length L}T_{, they form}

the total length L given by

L = LM + LT. (3.1)

In contrast to Hill’s original paper, muscle extension is defined as positive along the muscle length. In the case of an isometric contraction, in which the total length is kept constant, the muscle length LM _{will decrease while the tendon length L}T _{increases as the muscle is pulling on the tendon.}

The speed at which a muscle contracts is given by the rate of change of its total length with respect to time. From Eq. (3.1), the contraction speed is therefore

˙

L = ˙LM + ˙LT, (3.2)

where a superscribed dot indicates the time derivative.

3.1.2 Force Balance

From the free-body diagram of the Hill model shown in Fig.3.2, the force equilibrium leads to the following equations:

FT = FM+ FP E (3.3)

F = FT, (3.4)

(16)

CHAPTER 3. MUSCLE MODELS 9

where the externally applied force F equals the force in the tendon. It should be noted that for simplicity, it is common practice to consider the muscles to be massless and to add the muscle mass to the skeletal bone inertia matrix instead (Holmberg,2012). By this, we obtain static equilibrium for the muscle model, i.e. P_~

F = 0, leading to Eqs. (3.3) and (3.4).

FT

FM

FP E

F FT

Figure 3.2: Free-body diagram of the Hill model.

3.1.3 Muscle Force

The force generated by the muscle is not only dependent on the muscle length and the speed at which the muscle is contracted, it also depends on the amount of cross-bridges that are formed. While the length dependence takes into account the amount of active cross-bridges due to the filament sliding, the calcium dependence considers the amount of active cross-bridges based on the present calcium concentration, as more myosin heads can attach at the same time at a higher calcium concentration. This leads to an equation for the muscle force including a length dependence fL(LM), a speed dependence fv( ˙LM), and a dependence on the fraction of activated cross-bridges

φ:

FM = F₀M · fL LM · fv L˙M · φ. (3.5)

The equations for the length and speed dependences have to be chosen in a way that they satisfy experimental data, and several choices are possible.

The fraction of activated cross-bridges related to the calcium ions released into the sarcoplasm can be given by a first-order kinetic equation (Huxley, 1957). For simplicity, however, in the model, φ will be treated as a simple on/off switch that is either 1 (muscle activation) or 0 (muscle deactivation): φ = ( 0, deactivation 1, activation. (3.6)

3.1.4 Force-Length Dependence

The muscle force-length dependence differs from muscle to muscle. In general, it is bell-shaped and can sometimes be approximated by a Gaussian curve, see Fig. 3.3. The normalised length dependence fL(LM) in [0, 1] can therefore be approximated in a simple way by a Gaussian function:

fL LM = exp −

LM _{− µ}2 2σ2

!

, (3.7)

with µ being the mean and σ the standard variation.

In case of a force-length dependence as in Fig.2.4, this approximation does neither accurately represent the steep ascent for lengths close to the minimum length LM_min, nor the linear increase of force for larger lengths, as can be seen in Fig.3.3. However, LMopt as a distinctive parameter can

easily be fitted to different muscles.

Although the Gaussian approximation is easy to implement, it comes at a cost. The force-length dependence is not symmetric but has a longer tail to the right. In addition, a Gaussian distribution does never truly become zero, which is why a sensible ‘cut-off’ has to be chosen. Nonetheless, the advantage of fewer parameters, easy implementation, and differentiability makes the Gaussian distribution the preferred function herein.

(17)

FM

LM

F0M

LM

min LMopt LMmax

Experimental data Gaussian approximation

Figure 3.3: Gaussian force-length dependence compared to experimental data for the brachialis muscle presented inChang et al.(1999).

3.1.5 Force-Speed Dependence

The constant β in Eq. (2.1) can be calculated by setting the speed ˙LM _{to the maximum contraction}

speed ˙LM

max, which at the same time means that the muscle force FM is zero (see Fig. 2.6):

β = α · ˙L M max FM 0 . (3.8)

To write the Hill equation in a dimensionless form, the obtained β is then substituted into Eq. (2.1) and the Hill equation is solved for the muscle force FM normalised by the isometric force F0M:

fv= FM FM 0 = α FM 0 1 − ˙ LM ˙ LM max ˙ LM ˙ LM max + α FM 0 . (3.9)

As can be seen, Eq. (3.9) does only depend on the shortening speed ˙LM, as ˙LM_max and α/F₀M are constant. With k = α/F₀M, fv in [0, 1] can therefore also be seen as the speed dependence fv( ˙LM)

that is used to calculate the contractile force FM:

fv ˙LM = k 1 − ˙ LM ˙ LM max ˙ LM ˙ LM max + k . (3.10)

The maximum muscle speed ˙LM

max is determined by the type of fibres the muscle is composed of.

According to Nigg and Herzog (2006, p. 204), muscles composed of predominantly slow twitch fibres reach a maximum speed of ˙LM

max≈ 6 · LMoptm/s, while muscles composed of predominantly

fast twitch fibres are able to produce a maximum speed of ˙LM

max ≈ 16 · LMopt m/s. Most human

skeletal muscles are composed of both fibre types, with one of them being predominant, which is why the percentage of each fibre type has to be assessed.

3.1.6 Parallel Element and Tendon Force

Both the parallel and the serial element have an exponential force-deformation response over their respective lengths, LM _{and L}T_{. Several ways to model these forces exist, see e.g.}_Winters₍₁₉₉₀_).

3.2 Modified Skeletal Muscle Model

Setting up a Hill muscle model generally leads to more unknowns than equations, which is tradi-tionally solved by setting the tendon speed ˙LT _{to zero. The model presented in Fig.}_3.4 _{takes a}

(18)

different approach to muscle modelling by introducing an evolution law for LM_{. With the evolution}

law, LM _{and the contracting speed ˙}_LM _{can be numerically calculated considering the speeds of}

the muscle and tendon.

The contractile element is modelled by a friction clutch in parallel to a viscous damper. The friction clutch symbolically represents the cycling of the myosin heads as they attach to the thin filaments and pull them towards the centre of the sarcomere. The choice of using a friction clutch is motivated by the macroscopic motion of the filaments in contraction. Debold et al.(2005) showed that overall forward movement of the actin filament also included short and rapid backward slips, resulting in a microscopic jerking motion. Although some cross-bridges are always pulling on the thin filaments during contraction, the filaments will slide back slightly when myosin heads deattach during the contraction cycle, a behaviour which can be modelled by a friction clutch (Sharifimajd and St˚alhand,2013).

To account for a muscle’s viscoelasticity, a viscous damper is included parallel to the friction clutch and passive element. The viscosity is thought to arise from the fluid in the muscle resisting contraction, which is why the viscous damper can be thought of as part of the contractile element (McMahon, 1984, p. 14). Since connective tissue contains a lot of collagen, the parallel element is modelled with an exponential function. The corresponding free-body diagram of the modified model is shown in Fig.3.5.

F F PE LT LM M D T ω Mt CE L

Figure 3.4: Modified Hill model with friction clutch (M), damper (D), parallel element (PE), and tendon (T).

3.2.1 Kinematics

As can be seen in Fig.3.4, the total muscle length is given by the length of the contractile element, LM, and the tendon length LT,

L = LM + LT, (3.11)

with its time derivative being the contraction speed given by ˙

L = ˙LM + ˙LT. (3.12)

3.2.2 Force Balance

From the force equilibrium in the free-body diagram (Fig.3.5), and assuming static equilibrium as in section3.1.2, the tendon force FT is given by the sum of the muscle force FM, the parallel element force FP E, and the force of the viscous damper FD:

(19)

where the external force F of the model is equal to the tendon force

F = FT. (3.14)

The input of power used to generate muscle contraction is associated with the externally applied torque Mt:

Mt= rFM, (3.15)

where r is a constant. The power added through Mtω is thought to be the mechanical equivalence

to the power released by splitting ATP in the cross-bridge cycle.

3.2.3 Dissipation Inequality

The external power of the model is given by

P = F ˙L + Mtω, (3.16)

which, with Eqs. (3.12), (3.14) and (3.15), leads to

P = FT ˙LM + ˙LT+ FMrω. (3.17) Using Eq. (3.13) and rearranging the terms, we obtain

P = FTL˙T + FP EL˙M + FDL˙M + FM ˙LM + rω. (3.18) Disregarding all thermal effects, the free energy of the system becomes the strain-energy Ψ in the springs. It seems reasonable to assume that the total strain-energy is the sum of the strain-energy in the parallel and serial spring

Ψ LM, LT = ΨP E _LM_{+ Ψ}T _LT_. _(3.19)

The second law of thermodynamics requires that the change in elastic energy ˙Ψ cannot exceed the externally supplied power, i.e.

˙

Ψ ≤ P. (3.20)

The strain-energy is dependent on both the muscle length LM _{and the tendon length L}T _which

themselves are dependent on the time t. Using the chain rule, ˙ Ψ =dΨ L M_{, L}T dt = ∂Ψ ∂LML˙ M ₊ ∂Ψ ∂LTL˙ T_. _(3.21)

Substituting Eqs. (3.18) and (3.21) into (3.20), this leads to ∂Ψ ∂LML˙ M₊ ∂Ψ ∂LTL˙ T _{≤ F}T_L_˙T _{+ F}P E_L_˙M_{+ F}D_L_˙M _{+ F}M ˙LM _{+ rω}_, FT FM FP E F FT FD ω Mt

(20)

and after rearranging the terms FP E+ FD− ∂Ψ ∂LM ˙ LM+ FT − ∂Ψ ∂LT ˙ LT + FM ˙LM + rω≥ 0. (3.22) If the serial spring is assumed to be elastic, the tendon force can be expressed by

FT = ∂Ψ

∂LT. (3.23)

Back-substituting Eq. (3.23) and Eq. (3.13) into Eq. (3.22) gives ∂Ψ ∂LT − F M − ∂Ψ ∂LM ˙ LM + FM ˙LM + rω≥ 0. (3.24) This inequality must be satisfied for all muscle contractions.

3.2.4 Muscle Force

The tangential speed of the friction clutch is the product of its radius r and its angular velocity ω. The radius and the angular velocity of the friction clutch are unknown, their product can however be interpreted as the maximum contraction speed

˙

LM_max= rω. (3.25)

This interpretation is reasonable, as the actin cannot translate faster than the myosin heads cy-cle, and is therefore physiologically restricted by the muscle-specific maximum contraction speed

˙ LM

max. The muscle force FM is assumed to be proportional to the relative speed ( ˙LM+ ˙LMmax)

be-tween the friction clutch and the actin filament, which is portrayed in Fig.3.6. Remember that ˙LM

is positive in extension and therefore depicted in opposite direction of the contracting muscle force. In order for Eq. (3.24) to be satisfied, both addends are independently assumed to be greater or equal to zero. As shown in Fig.3.6, it is assumed that the force of the contractile element generated by the friction clutch is dependent on the relative speed ( ˙LM + ˙LMmax). Since no contractile force

is generated below the maximum contraction speed − ˙LMmax, the contractile force is set to zero for

˙

LM = − ˙LMmax. Thus, we assume

FM = (

κ ˙LM+ ˙LM_max, for ˙LM ≥ − ˙LM_max

0, otherwise,

(3.26)

where κ > 0, see Fig.3.7. The muscle force in Eq. (3.26) is linear in ˙LM _{but, to first order, resembles}

the speed dependence in Fig.2.6given by the Hill equation. This is a simplified assumption that

˙ LM ω ˙ LM max Myosin

Actin Binding sites

Figure 3.6: Illustration of the kinematics of the friction clutch model the cycling of the myosin heads (left image reproduced fromTortora and Derrickson,2013, p. 336).

(21)

can be adjusted if needed under the condition that Eq. (3.24) is satisfied. A better fit to the Hill curve is obtained by taking, e.g.,

FM =    κ ˙LM_{+ ˙}_LM max 3 , for ˙LM _{≥ − ˙}_LM max 0, otherwise. (3.27)

In order to account for the filament overlap, i.e. the length dependence fL(LM), and the calcium

dependence φ of the muscle force, the parameter κ is defined to be

κ = ˆκ · fL LM · φ, (3.28)

where ˆκ > 0 is a muscle-specific parameter. As described in section 3.1.4, fL(LM) is given by

Eq. (3.7), while the calcium dependence φ will be treated as a simple on/off switch according to Eq. (3.6). Since ˆκ has to apply to all FM_{, it is determined at the maximum isometric force F}M 0

at which ˙LM _{= 0, f}

L(LM) = 1, and φ = 1. Using Eqs. (3.26) and (3.28), follows

ˆ κ = F M 0 ˙ LM max . (3.29) − ˙LM max κ ˙LMmax ˙ LM FM 0

Figure 3.7: Muscle force FM _{from Eq. (}_3.26_{) over muscle speed ˙}_LM_.

3.2.5 Evolution of L

M

Following the same approach as in section3.2.4, the first addend in Eq. (3.24) is assumed to be FT − ∂Ψ

∂LM − F

M _{= ζ ˙}_LM_, _(3.30)

where ζ > 0. Substituting Eqs. (3.26) into (3.30) for ˙LM ≥ − ˙LM_max leads to ζ ˙LM + κ( ˙LM + ˙LM_max) = ∂Ψ

∂LT −

∂Ψ

∂LM, (3.31)

which rearranged gives an evolution law for LM ˙ LM = 1 ζ + κ ∂Ψ ∂LT − ∂Ψ ∂LM − κ ˙L M max . (3.32)

3.2.6 Specilisation to Hill’s Equation

Equation (3.31) can be further specialised to recover the Hill equation in Eq. (2.1). Towards this end, take the terms on the left side of Eq. (3.31) to be

(22)

Using Eq. (3.31) together with Eq.(3.23), the afterload τ can be defined to be τ = ∂Ψ ∂LT − ∂Ψ ∂LM = F T₋ ∂Ψ ∂LM. (3.34)

For an isometric contraction, the length LM _{is constant at steady state, which is why ˙}_LM _{= 0 and}

therefore, applied to Eq. (3.33),

τ0= κ ˙LMmax. (3.35)

The stress τ0can be interpreted as the maximum contractile force depending on the current muscle

length LM_{. Eq. (}_3.33_{) can thus be rewritten to}

˙

LM(ζ + κ) = τ − τ0. (3.36)

In order to end up with the Hill equation (Eq.2.1), define ζ + κ = 1

β(τ + α), (3.37)

where α, β > 0. With Eq. (3.36) this finally leads to ˙

LM = βτ − τ0

τ + α. (3.38)

It should be noted that contraction is defined negatively and extension positively in Eq. (3.38), which is why the term τ − τ0 has opposite signs from the original Hill equation (Eq.2.1).

According toHill(1938), the ratio α/FM

0 = 0.25 and from Eq. (3.8) follows:

α = 0.25 · F₀M (3.39)

(23)

4 Computer Implementation

In a first step, two quasi-static simulations modelling the change in length during muscle con-traction are implemented in Matlab version 9.4.0.813654 (R2018a), using the equations stated in the previous chapters. The quasi-static simulations are performed under isometric and isotonic conditions, respectively. They are based on the modified Hill model described in section3.2and simulate the change in total length, taking into account that the tendon length is not constant. Following the experimental protocol inB¨ol et al.(2013), these simulations are performed to verify the introduced modified muscle model.

The quasi-static simulations will subsequently be merged into a more accurate simulation that is neither isometric nor isotonic, as these conditions are not realistic during actual muscle contraction and body movement. This general contraction simulation models a biceps curl for the brachialis muscle for a given path of motion and speed with no additional weight. The same biceps curl is also modelled in AnyBody, which assumes a constant tendon length. The results of both biceps curl simulations will eventually be compared to each other in order to study the magnitude of the errors that are introduced by assuming a constant tendon length.

4.1 Physiological Muscle Parameters

Each muscle has a set of distinctive physiological parameters that have to be considered in the simulation. These parameters vary from specimen to specimen and are sometimes taken from ani-mals and extrapolated to humans. The parameters for the brachialis muscle used in the AnyBody and in the Matlab simulations are collected in Table4.1.

The physiological cross-sectional area (PCSA) is given by the ratio of the muscle volume VM

to the muscle belly length LM

B, corrected for the pennation angle γ, and is proportional to the

maximum isometric force FM

0 (Zajac,1989): P CSA = V M LM B cos(γ) . (4.1)

The corresponding tendon strain due to maximum isometric force is given by εT

0. The tendon slack

length LT_s determines the tendon length at which the tendon starts producing force. The column ‘Fibre type’ indicates whether a muscle is composed of predominantly slow twitch (ST) or fast twitch (FT) muscle fibres. Knowing the percentage of each fibre type and together with the values given in section3.1.5, the maximum contraction speed ˙LMmax is defined to be

˙ LM_max= 6 · LM_opt·%ST 100s + 16 · L M opt· %F T 100s. (4.2)

According toGarner and Pandy (2003), the minimum and maximum muscle lengths LM min and

LM

maxare taken to be

LM_min = 0.5 · LM_opt (4.3)

LMmax= 1.5 · LMopt. (4.4)

It should be noted that for simplicity, the pennation angle is set to zero in the AnyBody as well as the Matlab simulations, despite its value being γ = 15◦ _{for the brachialis muscle according}

to the AnyBody Managed Model Repository.

(24)

CHAPTER 4. COMPUTER IMPLEMENTATION 17 VM_[ml] _LM B [m] LMopt[m] γ [◦] LTs [m] εT0 [−] JT[−] Fibre type 0.5 · 102 0.173 0.123 0 0.15 · LMopt 0.053 3.0 40% FT

Table 4.1: Muscle parameters for the brachialis muscle taken from the AnyBody Managed Model Repository v. 1.6.3 (γ = 15◦ originally).

4.2 Parallel Element and Tendon Force

As mentioned in section3.1.6, the parallel element and tendon force can be modelled with expo-nential functions suited for the muscle in question. In the following Matlab simulations, the tendon force is chosen to be specified as FT directly rather than the strain-energy in the serial spring ΨT. As opposed to what is shown in Fig. 2.5, where the tendon force is dependent on the muscle length LM_{, the tendon force exhibits an exponential curve when plotted over the tendon length L}T

(Zajac, 1989). Therefore, during contraction (LT _{≥ L}T

s), the tendon force can be approximated

by an exponential function that is defined by FT _{= 0 at L}T _{= L}T

s, and FT = F0M at LT = LT0

(Christensen et al.,2002). The tendon length at which the tendon force equals the muscle-specific maximum isometric force is given by

LT₀ = LT_s 1 + εT₀ . (4.5)

The tendon force is modelled according to the following equation proposed byChristensen et al.

(2002): FT = F M 0 exp (JT_{) − 1} exp JTL T _{− L}T s LT 0 − LTs − 1 , (4.6)

with JT _{as a dimensionless shape parameter whose value is given in Table}_4.1_{. The strain-energy}

function of the parallel spring is taken to be (Holzapfel et al.,2000) ΨP E = C1 2C2 exp C2 λM 2− 12 − 1 , (4.7)

where C1and C2are material parameters specific to the muscle that is being modelled, and λM ≥ 1

is the stretch of the muscle given by the ratio of the current muscle length LM _{to its optimal length}

LM opt, i.e. λM = L M LM opt . (4.8)

The derivative dΨP E/dLM can be calculated with Eqs. (4.7) and (4.8), giving dΨP E dLM = 2C1 λM LM opt λM 2− 1exp C2 λM 2− 1 2 . (4.9)

4.3 Quasi-Static Simulations

Quasi-static conditions are assumed herein. The chosen starting point for the quasi-static con-traction simulations is LM

opt, which is the muscle length at which the parallel element force FP E is

assumed to be zero (see Fig.2.5). By choosing LM

opt, the stretch of the muscle λM will become 1

and dΨP E_/dLM _{therefore zero (Eqs.} _4.8 _and_4.9_{). Additionally, it is the muscle length at which}

the muscle produces its maximum force and is therefore a value of interest during isometric and isotonic muscle experiments.

(25)

CHAPTER 4. COMPUTER IMPLEMENTATION 18

4.3.1 Isometric Contraction

During an isometric contraction, the total length L is kept constant, i.e. the muscle length LM

decreases while the tendon length LT _{increases in such a way that the total length is maintained.}

The stress τ0 represents the maximum isometric muscle force, qualified by the length

depen-dence fL(LM), and changes during contraction. The length dependence of the brachialis muscle is

modelled by a Gaussian function (Eq.3.7) which approximates the reality reasonably well (Chang et al., 1999). In Matlab, the Gaussian function is fitted to the optimal, minimum and maximum muscle length by using the non-linear least-squares solver lsqcurvefit. The found parameters p(1) and p(2) represent the mean µ and the standard deviation σ, respectively, by which the Gaussian length dependence is defined. In the case of a simple Gaussian function, the mean can be manually set to the optimal length which allows for the standard deviation to be calculated analytically. However, in the current state, the code can be easily adjusted to fit a more complex function with additional parameters, e.g. a skewed Gaussian function. Eq. (3.7) can be further adjusted to match the length dependence given byChang et al.(1999), see Fig.3.3.

1 x d a t a = [ L M m i n L M o p t L M m a x ]; 2 y d a t a = [0 1 0]; 3 4 % p (1) e q u a l s mean , p (2) e q u a l s s i g m a 5 G a u s s = @( p , x d a t a ) exp ( -( xdata - p (1) ) . ^ 2 / ( 2 * p (2) ^2) ) ; 6 p = l s q c u r v e f i t ( Gauss ,[ L M o p t 0.1] , xdata , y d a t a ) ; 7 f L _ x d a t a = l i n s p a c e ( LMmin , L M m a x ) ; 8 fL = G a u s s ( p , f L _ x d a t a ) ;

In order to calculate the muscle length LM_{, the Matlab ordinary differential equation (ODE)}

solver ode45 is used, which determines the solution iteratively over a given period of time and with a given initial condition. All used equations in the ODE solver are made dependent on LM _as

the only unknown variable, using Eq. (3.11). Following the experimental protocol for an isometric contraction inB¨ol et al.(2013), the muscle length LM _{is brought to the optimal length L}M

optbefore

activation. It is therefore the initial muscle length fed to the ODE solver. The Matlab commands sol and deval are used to evaluate the solution of ode45 at any point in time.

1 d L M d t = @( t , LM ) b e t a *( tau ( LM ) - t a u 0 ( LM ) ) /( tau ( LM ) + a l p h a ) ;

2 sol = o d e 4 5 ( dLMdt , time , L M o p t ) ;

3 [ LM , L M d o t ] = d e v a l ( sol , t i m e ) ;

Activation and deactivation of the muscle are driven by the calcium dependence activation switch phi (Eq.3.6), only affecting tau0. The stresses tau and tau0 are computed from Eqs. (3.34) and (3.35), disregarding the strain-energy which is zero. The constants alpha and beta are given in Eqs. (3.39) and (3.40), respectively. The resulting muscle length is then subtracted from the constant total length to obtain the tendon length LT_{, which is used in Eq. (}_4.6_{) to calculate the}

outgoing force of the model, the tendon force FT_.

An overview of the relations of the equations to one another and their implementation in Matlab is given in Fig.4.1.

4.3.2 Isotonic Contraction

During an isotonic contraction, the externally applied force of the model is kept constant. In a real experiment this is often achieved by letting the muscle contract against a hanging weight. The external force is therefore

F = W, (4.10)

with W being the weight. Following Eq. (3.4), the force F is at the same time the tendon force FT_{, leading to a constant tendon length L}T _{during an isotonic contraction if inertial effects are}

(26)

neglected, as the constant tendon force only depends on the tendon length. The muscle length LM

and therefore also the total length L will however change during the activation of the muscle. At a constant tendon force and no contraction of the muscle, the tendon length will assume an isotonic tendon length LT

isot, which will not change during isotonic activation, since quasi-static

conditions are assumed. The isotonic tendon length is determined via non-linear root finding by setting the tendon force in Eq. (4.6) equal to the outer force (Eq. 4.10) and solving for LT. In Matlab, this is accomplished by using the command fzero, which will find a tendon length at which the error is minimised. The tendon slack length LTs is used as an initial value to the function.

1 e r r o r = @( LT ) F - F0 /( exp ( JT ) -1) *( exp ( JT .*( LT - LTs ) /( LT0 - LTs ) ) -1) ;

2 L T _ i s o t = f z e r o ( error , LTs ) ;

Following the experimental protocol for an isotonic contraction inB¨ol et al.(2013), the muscle length is first brought to the optimal muscle length LM

opt while keeping the tendon length at the

isotonic tendon length LT

isot. It should be noted that in the actual muscle experiments conducted

byB¨ol et al.(2013), the change in muscle length provokes a small tendon force before the tendon force is kept constant during activation, as the muscle is stretched from its resting length to the optimal muscle length and will therefore pull on the tendon. In this simulation, however, the tendon length will be assumed to stay at its isotonic tendon length while the muscle length starts at its optimal length, not provoking a tendon force. The first change in tendon force is provoked when the muscle is activated.

The Gaussian force-length dependence fL(LM) of the muscle force is fitted to the physiological

data as done during the isometric contraction.

Despite the isotonic contraction starting at the optimal muscle length, the strain-energy (Eq.4.9) simulating the passive force will be required during deactivation, as opposed to the isometric con-traction. With no muscle force present, the muscle will elongate due to the attached weight until the stretch of the muscle induces a passive force great enough to withstand elongation. Similarily to the Gaussian function, the parameters of the strain-energy function are obtained using the non-linear least-squares solver lsqcurvefit. The strain-energy is fitted to the parallel element force Fpe over the muscle length Lmobtained from an AnyBody simulation.

1 x d a t a = L m _ a n y b o d y ; 2 y d a t a = F p e _ a n y b o d y ; 3 4 S t r a i n _ e n e r g y = @( C , x d a t a ) 2* C (1) *( x d a t a /( L M o p t ^2) ) . * ( ( x d a t a / L M o p t ) .^2 -1) .* exp ( C (2) . * ( ( x d a t a / L M o p t ) .^2 -1) . ^ 2 ) ; 5 C = l s q c u r v e f i t ( S t r a i n _ e n e r g y ,[1 1] , xdata , y d a t a ) ; Since the muscle is brought to the optimal muscle length LM

optand subsequently activated, LMopt

is used as the initial condition for calculating the muscle length. As the total length is not given in contrast to isometric contractions, both LM _{and L are unknown during the computation of L}M_.

The muscle length and its time derivative during activation and deactivation are then calculated as done for the isometric contraction using the ODE solver ode45, where tau0 is set to zero through the calcium concentration phi during deactivation of the muscle.

Using the obtained muscle speed ˙LM_{, the muscle force is eventually calculated by Eq.}_3.26_.

4.4 General Contraction Simulations

A more realistic simulation that is neither isometric nor isotonic is performed in Matlab, as both are simplified assumptions that do not hold during actual muscle movements. Quasi-static conditions will however be assumed for the muscle force FM and the parallel element force FP E, but not

(27)

for the tendon force FT_{, as bones have inertia. In the following, parameters computed in Matlab}

are denoted by superscribed capital letters, while lengths and forces belonging to AnyBody are denoted by subscribed lowercase letters.

In AnyBody, the muscle length Lm is calculated by setting the tendon speed ˙Ltto zero. The

simulation for a simple biceps curl is repeated in Matlab using the same physiological parameters and some of the values obtained from the AnyBody simulation to compute the muscle length LM without setting the tendon speed ˙LT to zero. Among others, the calculated muscle forces FM and Fm are subsequently compared to each other.

To ensure comparability, the Matlab simulations will always have the same time duration as the AnyBody simulations. The input data to both the AnyBody and Matlab simulations are the total length Lmt, the tendon force Ft, and time. The parameters Lmt and ˙Lmt are taken to be

given, as the movement in the simulation is predetermined. In AnyBody, the tendon force Ft is

then calculated by the Newton-Euler equations.

The muscle force can be scaled with an activation variable, a, that determines the amount of voluntary contraction. It is given by the ratio of the muscle force Fmto the available muscle force,

the so-called Strength:

a = Fm

Strength. (4.11)

As the passive parallel element force will always be at its maximum and does not depend on the activation, it is energetically favourable for the muscle to decrease its activation if the passive force will take over part of the external load. A simulation in AnyBody will therefore always try to minimise the activation a if a passive force is present. The equivalent of the activation a is the calcium dependence φ in Matlab, which is assumed to be 1 at all times. As AnyBody provides the activation a and also the Strength in their results, the fully activated muscle force FM _{in Matlab}

is eventually compared to the Strength in AnyBody.

In order to approximate the force-length dependence used in the AnyBody simulation, the Gaussian distribution (Eq.3.7) is fitted to the Strength obtained from AnyBody normalised by the maximum isometric force FM

0 using the non-linear least-squares solver lsqcurvefit.

1 S t r e n g t h _ a n y _ n o r m = S t r e n g t h _ a n y b o d y / F_0 ; 2 3 % p (1) e q u a l s mean , p (2) e q u a l s s i g m a 4 G a u s s = @( p , L m _ a n y b o d y ) exp ( -( L m _ a n y b o d y - p (1) ) . ^ 2 / ( 2 * p (2) ^2) ) ; 5 p = l s q c u r v e f i t ( Gauss ,[ L M o p t 0.01] , L m _ a n y b o d y , S t r e n g t h _ a n y _ n o r m ) ; 6 x D a t a = t r a n s p o s e ( l i n s p a c e ( LMmin , LMmax , l e n g t h ( t i m e _ s i m ) ) ) ; 7 fL = G a u s s ( p , x D a t a ) ;

The same approach was chosen for the strain-energy of the parallel spring dΨP E_/dLM _(Eq._4.9_),

which is fitted to the parallel element force Fpeto obtain the parameters C1and C2, identically to

the isotonic contraction simulation.

1 x d a t a = L m _ a n y b o d y ; 2 y d a t a = F p e _ a n y b o d y ; 3 4 S t r a i n _ e n e r g y = @( C , x d a t a ) 2* C (1) *( x d a t a /( L M o p t ^2) ) . * ( ( x d a t a / L M o p t ) .^2 -1) .* exp ( C (2) . * ( ( x d a t a / L M o p t ) .^2 -1) . ^ 2 ) ; 5 C = l s q c u r v e f i t ( S t r a i n _ e n e r g y ,[1 1] , xdata , y d a t a ) ;

As done in the isometric and isotonic contraction simulations, the muscle length LM _{and the}

muscle speed ˙LM _{will be calculated using the ODE solver ode45. All functions used in the ODE}

solver are made dependent on LM _{and L}

mt, using Eq. (3.11). The first value for LM that is fed

(28)

externally applied tendon force Ft at the first time step, and solving for LT. In Matlab, this is

accomplished by using the non-linear root finding function fzero.

1 LT0 = LTs * ( 1 + s t r a i n _ T 0 ) ; 2 3 e r r o r = @( LT ) F t _ a n y b o d y (1) - F0 /( exp ( JT ) -1) *( exp ( JT *( LT - LTs ) /( LT0 - LTs ) ) -1) ; 4 L T _ i n i t i a l = f z e r o ( error , LTs ) ; 5 LM (1) = L m t _ a n y b o d y (1) - L T _ i n i t i a l ; To obtain one LM _{to each given L}

mt, the ODE function is looped through the vector length of

the time of the simulation, time sim. However, since the muscle speed ˙LM _{is calculated based on}

two values of the muscle length LM_{, the vector of L}M _{needs to be greater than the vector of ˙}_LM

by one value in order to adequately calculate ˙LM_{. The time of the simulation used in the ODE}

solver is therefore artificially prolonged in a way consistent with the time steps in time sim. The maximum and last value of the vector time sim is the length of the simulation in seconds. The last entry of the muscle length LM is deleted after the ODE solver is exited in order to obtain equally long vectors.

1 t i m e _ o d e = t i m e _ s i m ; 2 t i m e _ o d e ( end +1) = t i m e _ s i m ( end ) + t i m e _ s i m (2) ; 3 4 for i = 1:( l e n g t h ( t i m e _ o d e ) ) 5 6 FT = @( LM ) F0 ./( exp ( JT ) -1) .*( exp ( JT .*( L m t _ a n y b o d y ( i ) - LM - LTs ) ./( LT0 - LTs ) ) -1) ; 7 fL = @( LM ) exp ( -( LM - p (1) ) . ^ 2 / ( 2 * p (2) ^2) ) ; 8 t a u 0 = @( LM ) k a p p a _ h a t * v _ m a x * fL ( LM ) * phi ; 9 d p s i _ d L M = @( LM ) 2* C (1) *( LM /( L M o p t ^2) ) *(( LM / L M o p t ) ^2 -1) * exp ( C (2) *(( LM / L M o p t ) ^2 -1) ^2) .*( LM / LMopt > = 1 ) ; 10 tau = @( LM ) FT ( LM ) - d p s i _ d L M ( LM ) ; 11 12 t i m e _ a c t = [ t i m e _ o d e ( i ) , t i m e _ o d e ( i +1) ]; 13 d L M d t _ a c t = @( t , LM ) b e t a *( tau ( LM ) - t a u 0 ( LM ) ) /( tau ( LM ) + a l p h a ) ; 14 sol = o d e 4 5 ( d L M d t _ a c t , t i m e _ a c t , LM ( i ) ) ; 15 [ LM_act , L M d o t _ a c t ] = d e v a l ( sol , t i m e _ a c t ) ; 16 LM ( i +1 ,1) = L M _ a c t (2) ; 17 L M d o t ( i ) = L M d o t _ a c t (1) ; 18 L M d o t ( i +1) = L M d o t _ a c t (2) ; 19 20 if l e n g t h ( LM ) == l e n g t h ( t i m e _ o d e ) 21 b r e a k 22 end 23 end 24 25 LM ( end ) = [ ] ;

The tendon length LT and the tendon speed ˙LT are calculated using the kinematic relations in Eqs. (3.11) and (3.12), respectively. The available muscle force FM can subsequently be calculated using Eq. (3.26).

(29)

4.4.1 Modified Coefficient of Determination

To be able to compare the curves of the lengths, speeds and forces in AnyBody and Matlab, a modified coefficient of determination ˆR2 is introduced. Instead of calculating the established coefficient of determination R2, the mean of the respective curves is subtracted before calculating R2_{. By correcting the curves by their mean difference, the shape of the curves is compared rather}

than their actual values.

4.4.2 Muscle Work

The area under the muscle force curves is equivalent to the muscle work and is calculated using the Matlab function trapz. The muscle work provides the energy used to generate the muscle force and is of physiological interest.

1 A r e a _ F M = abs ( t r a p z ( LM , FM ) ) ;

4.4.3 Error Percentage

In order to estimate the magnitude of the errors introduced by assuming a constant tendon length, the deviation of the calculated tendon speed to the total speed is regarded. This is done by calculating the ratio of ˙LT _{to ˙}_{L. This ratio will always be zero in AnyBody, as the tendon speed}

˙

Lt is zero at all times. While the total length and therefore the total speed is identical in both

AnyBody and Matlab, the split between the muscle and tendon according to Eq. (3.11) will differ in Matlab, and therefore also the ratio of ˙LT to ˙L.

(30)

CHAPTER 4. COMPUTER IMPLEMENTATION 23 τ0 = κ ˙ L M max ˆκ = F M 0 ˙ L M max φ F M 0 ε T 0 L T s ˙ L M max ˙ L M (i ) ˙ L M L M L M (i+1 ) ode45 End? Exit L T 0 = L T s 1 + ε T 0 κ = ˆκ · fL L M · φ fL L M = exp − L M − µ 2 2 σ 2 ! τ = F T − ∂ Ψ ∂ L M ˙ L M = β τ − τ0 τ + α F T = F M 0 exp (J T) − 1 exp J T L T − L T s L T 0 − L T s − 1 Input P arameters from literature Y N i=i+1 L M (i ) F T = F M 0 exp (J T) − 1 exp J T L T − L T s L T 0 − L T s − 1 L M opt L L T = L − L M J T α = 0 .25 · F M 0 β = 0 .25 · ˙ L M max

(31)

CHAPTER 4. COMPUTER IMPLEMENTATION 24 L M opt τ0 = κ ˙ L M max ˆκ = F M 0 ˙ L M max φ F M 0 ε T 0 L T s J T W ˙ L M max L M (i ) ˙ L M (i ) ˙ L M L M F M = κ ˙ L M + ˙ L M max α = 0 .25 · F M 0 β = 0 .25 · ˙ L M max L M (i+1 ) fzero ode45 Y End? Exit L T 0 = L T s 1 + ε T 0 κ = ˆκ · fL L M · φ fL L M = exp − L M − µ 2 2 σ 2 ! d Ψ P E d L M = 2 C1 λ M L M 0 λ M 2 − 1 exp C2 λ M 2 − 1 2 λ M = L M L M opt τ = F T − ∂ Ψ ∂ L M ˙ L M = β τ − τ0 τ + α F T = F M 0 exp (J T) − 1 exp J T L T − L T s L T 0 − L T s − 1 N i=i+1 Fpe C1 , C2 P arameters from literature L T isot Input Output from An yBo d y P arameters from literature

(32)

CHAPTER 4. COMPUTER IMPLEMENTATION 25 L M opt τ0 = κ ˙ L M max ˆκ = F M 0 ˙ L M max Fpe φ Fm F M 0 ε T 0 L T s J T L ˙ L M max C1 , C2 L T = L − L M L M (i ) F T (t = t1 ) ˙ L M (i ) ˙ L M L M F M = κ ˙ L M + ˙ L M max α = 0 .25 · F M 0 β = 0 .25 · ˙ L M max L M (i+1 ) fzero ode45 Y End? Exit L T 0 = L T s 1 + ε T 0 κ = ˆκ · fL L M · φ fL L M = exp − L M − µ 2 2 σ 2 ! d Ψ P E d L M = 2 C1 λ M L M 0 λ M 2 − 1 exp C2 λ M 2 − 1 2 λ M = L M L M opt τ = F T − ∂ Ψ ∂ L M ˙ L M = β τ − τ0 τ + α F T = F M 0 exp (J T) − 1 exp J T L T − L T s L T 0 − L T s − 1 N Input Output from An yBo dy P arameters from literature L M (1 ) i=i+1

(33)

5 Results

In the following chapter, the results of the implemented friction clutch model are presented and compared to the results byB¨ol et al.(2013) and those obtained from AnyBody.

5.1 Quasi-Static Simulations

An isometric and an isotonic simulation were carried out to validate the introduced friction clutch model. The simulations followed the experimental protocol by B¨ol et al. (2013) and obtained similar force-length curves. These results are only meant to verify the shape of the curves and not necessarily their amplitudes, as different input data was used.

5.1.1 Isometric Contraction

As the muscle length LM and the tendon length LT are limited by a constant total length, an increase of the former will result in a decrease of the latter, and vice versa. During activation of the muscle, LM will decrease while LT increases, as shown in Fig. 5.1. During deactivation, the muscle length elongates while the tendon length must shorten at the same time in order to maintain the constant total length. The tendon force will therefore decrease with decreasing tendon length during deactivation, leading to a first-order length response when calculating the muscle length, as can be seen in Fig.5.2.

In contrast to Table4.1, the tendon slack length LT

s is taken to be 0.5·LMoptto avoid ripples in the

graphs. Moreover, the time intervals used during the isometric simulation are significantly shorter than the ones presented inB¨ol et al.(2013) in order to obtain similar results. The disparities might be related to the different muscles and the different parameters that have been used during the experiment and the simulation. In the Matlab simulation, the muscle is activated at t = 0.1 s and deactivated at t = 0.13 s. For the sake of completeness, it should be noted thatB¨ol et al.(2013) provides an additional step of bringing the muscle from its resting length to its optimal length before activation. This step has been disregarded, as it was deemed unnecessary for the purpose of the simulation. 0 0.05 0.1 0.15 0.2 0.25 0.3 0.05 0.06 0.07 0.08 0.09 0.1 0.11 0.12 0.13

Figure 5.1: Muscle length LM _and

ten-don length LT _{over time during an}

iso-metric contraction. 0 0.05 0.1 0.15 0.2 0.25 0.3 0 50 100 150 200 250 300

Figure 5.2: Tendon force FT _{over time}

during an isometric contraction.

(34)

CHAPTER 5. RESULTS 27

5.1.2 Isotonic Contraction

Following the reasoning in the previous section, the additional step of bringing the muscle from its resting length to its optimal length before activation as performed byB¨ol et al. (2013) has been disregarded in the isotonic contraction simulation. The simulation starts at LM = LMopt directly.

The muscle is activated at t = 1 s and, followingB¨ol et al.(2013), has an activation period of 1.2 s. The weight attached to the end of the muscle is set to W = 100 N.

During activation, the muscle length first decreases until it assumes a length corresponding to a muscle force that equals the external weight, as shown in Fig.5.3. Additional shortening of the muscle would lead to a muscle force decrease due to its force-length dependence, which in turn would lead to an elongation of the muscle due to the attached weight, resulting in an increased muscle force once again. The muscle length will therefore assume a value at which the muscle force and tendon force are in equilibrium. The peak in muscle force in the beginning of the activation period in Fig.5.4occurs due to the muscle being activated at its optimal muscle length, resulting in the maximum generated force. The muscle force decreases with further muscle shortening until it equals the tendon force of 100 N.

Since the total length L changes with the muscle length LM, the muscle length is not restricted by the kinematics during deactivation. The passive force restrains elongation for LM > LMopt

(Fig.5.4) and gives a first-order response as seen in Fig.5.3. As the passive force prevents further elongation of the muscle, it will assume the same value as the tendon force to maintain equilibrium, identical to the muscle force during activation.

0 0.5 1 1.5 2 2.5 3 3.5 4 0.1 0.12 0.14 0.16 0.18 0.2 0.22 0.24

Figure 5.3: Muscle length LM and to-tal length L over time during an isotonic contraction. 0 0.5 1 1.5 2 2.5 3 3.5 4 0 20 40 60 80 100 120 140 160 180 200 220

Figure 5.4: Muscle force FM and parallel element force FP E over time during an isotonic contraction.

5.2 General Contraction Simulations

Several simulations were conducted testing the friction clutch model under conditions that are neither isometric nor isotonic. Two of those are separate simulations belonging to one single contraction; one during the ascent and one during the descent of the muscle force curve. Two further simulations were performed during general shortening, i.e. general contraction, and one simulation during general lengthening, i.e. general elongation.

All calculated forces in Matlab and AnyBody are the maximum available forces, i.e. fully activated. While the AnyBody muscle force at a = 1 is denoted by Strength, FM will continue to denote the fully activated muscle force at φ = 1 in Matlab.

(35)

CHAPTER 5. RESULTS 28

5.2.1 Simulations During Ascent and Descent of Muscle Force Curve

The data obtained from AnyBody comes from two separate simulations of the same body move-ment, each being 1 s long. During the first simulation, the body movement was modelled during the ascent of the muscle force curve, i.e. LM < LMopt, while the second simulation modelled the descent

of the curve, i.e. LM _{> L}M

opt. Although the simulations were separated into ascent and descent,

the Gaussian force-length dependence fL(LM) was fitted simultaneously to both AnyBody output

Strengths, giving a relatively good fit that is, however, better on the ascent than on the descent part of the curve (Fig. 5.5). As mentioned in section4.4, the Gaussian distribution is fitted to AnyBody’s Strength normalised by the maximum isometric force FM

0 .

The calculated Matlab muscle force FM _{is around the same value and shape as the AnyBody}

Strength (Fig.5.6 to Fig.5.9). However, the curves for descent are less accurate because of the Gaussian force-length dependence having a better fit on the ascending slope.

The shape of the muscle length curves in Matlab and AnyBody are very similar (Fig.5.10and Fig.5.11), resulting in a modified coefficient of determination of ˆR2_{= 1.0000 for both ascent and}

0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0.22 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Figure 5.5: Gaussian force-length dependence fitted to the normalised ascent and descent of the Strength obtained in AnyBody.

0.094 0.096 0.098 0.1 0.102 0.104 0.106 0.108 0.11 0.112 150 160 170 180 190 200 210 220 230 240 250 260

Figure 5.6: Comparison of muscle forces in AnyBody and Matlab over their re-spective muscle length during ascent.

0.143 0.144 0.145 0.146 0.147 0.148 0.149 0.15 150 160 170 180 190 200 210 220 230 240 250 260

Figure 5.7: Comparison of muscle forces in AnyBody and Matlab over their re-spective muscle length during descent.

An Analysis of Including the Evolution Law for the Serial Element in the Musculoskeletal Modelling

Biomedical Engineering

Master’s Thesis