Muscles are responsible for generating the forces required for the movement of multicellular organisms. Microscopically, these forces arise as a consequence of motor proteins (myosin) pulling and sliding along actin filaments. Current knowledge states that the molecular forces between actin and myosin are linear in nature [Huxley and Simmons (1971) Nature (London) 233, 533–538] and that the physiologically observed non-linearities (e.g. Hill's force–velocity relationship) are a consequence of non-linearities in the attachment/detachment ratios. However, this view has been disputed recently [Nielsen (2002) J. Theor. Biol. 219, 99–119], inspired by results from protein pulling experiments showing that proteins often have non-linear entropic force–extension profiles. Irrespective of the case, the present study aims at integrating such basic force-producing properties into large-scale simulations of muscle, which may accommodate macroscopic properties of muscles, e.g. the catch-like effect, the Henneman principle and accurate twitch force and motor unit size distributions. As a test of the underlying principles, a model of the biceps caput breve muscle is presented and compared with experimental data.

Introduction

The complex relationship that exists between molecular muscle structure and function has generated a great need for accurate models of muscles. In the present study, an attempt is made to integrate a number of old and new biophysical models of muscle function, resulting in a model of whole muscle that presents many of the most important physiological characteristics while still maintaining a strict and close parametric relationship with the underlying hierarchies of structures all the way down to the molecular level.

Models and methods

The view that the tensile stress between myosin and actin may best be described using linear spring forces and non-linear attachment/detachment ratios as proposed by Huxley and Simmons [1] was opposed recently by Nielsen [2,3] on the basis of recently published force–extension profiles derived from protein pulling experiments [4,5]. Such experiments revealed ‘entropic elasticity’ as the underlying tension-generating mechanism. This notion was used by Nielsen [2] to find the velocity-dependent sarcomeric force Fsarc by deriving the following set of equations:

 
formula
(1)

where L is the contour length of the force-producing region of the myosin chain, kB is Boltzmann's constant, T the temperature and A the ‘persistence length’. The parameters a, b and N are linearly dependent on contraction velocity, where [a b] represents the interval of allowed cross-bridge lengths at a given velocity and N is the number of effective cross-bridges formed at a given velocity and for a given muscular activation level (see [2] for details).

In [6], it is assumed that the force-producing region of myosin may uncoil at high loads, causing changes in the contour length L, thus effectively creating two populations of myosin within the sarcomere with slightly different force-producing capabilities. According to Nielsen [6], a linear combination of these may account for the double hyperbolic force–velocity relation observed by Edman [7].

Regardless of exactly how the sarcomeric tension is calculated (e.g. by following Hill [8], Huxley and Simmons [1] or Nielsen [2]), the sarcomeric tension has to be scaled up to get the total force of a muscle. In the remainder of the present study, therefore, we shall be concerned mainly with the scaling factors that link the molecular forces responsible for the sliding movement of myosin and actin filaments with the forces observed in whole muscles.

Model of the BISH (biceps caput breve)

For the purpose of illustration, the following model is based on data from the BISH. As neural input, 800 regular spiking neurons were modelled according to the prescriptions given by Wilson [9] for regular spiking neurons. As suggested by Nielsen [3], a neuron-size-dependent modifier ρi was added that affects the input currents given to the individual neurons i by scaling the currents according to the surface area of the neurons. It has the form

 
formula
(2)

where carea is an arbitrary scaling constant and ri is the radius of the ith neuron. The neuron radii were distributed according to

 
formula
(3)

where c1 and c2 are arbitrary constants and i is the index to the neurons. In the present study, these constants were chosen such that ρi had values in the interval [0.08 0.8]. These neurons represent the α-motoneurons and were connected to an estimated nmf=2.7×108 myofibrils in BISH (assuming cross-sectional areas 2.1 cm2 and 7.8×105 nm2 for BISH and a sarcomere respectively).

The fraction φi of myofibrils assigned to the ith motor unit is

 
formula
(4)

where φnorm is a normalization factor, ∑i=1nμ·φi=1 and βμ depends linearly on the number of units in the muscle nμ. The resulting twitch force distribution in the motor pool can be seen in Figure 1(A).

Simulation results

Figure 1
Simulation results

(A) Motor unit composition in BISH muscle simulation. This exponential distribution of motor unit sizes corresponds very well to the experimental data reported by Milner-Brown et al. [11] and Monster and Chan [12], who concluded that an exponential relation exists linking the number of motor units to the motor unit forces. (B) Firing frequency of individual motor units pertaining to the BISH muscle simulation. These simulation results are found to be qualitatively quite similar to the experimental data reported in [12]. Note that each line corresponds to a motor neuron. (C) Twitch force of recruited fibres as a function of the total force produced in BISH muscle simulation. There is a good correspondence between this result and the data reported in [11,12]. (D) The total muscle force as a function of the percentage of motor units that have been recruited (in agreement with Walmsley et al. [13]). (E) Gradual increase in BISH muscle force during a ramped increase in input current to the monitor unit neurons. Modified from Nielsen [3].

Figure 1
Simulation results

(A) Motor unit composition in BISH muscle simulation. This exponential distribution of motor unit sizes corresponds very well to the experimental data reported by Milner-Brown et al. [11] and Monster and Chan [12], who concluded that an exponential relation exists linking the number of motor units to the motor unit forces. (B) Firing frequency of individual motor units pertaining to the BISH muscle simulation. These simulation results are found to be qualitatively quite similar to the experimental data reported in [12]. Note that each line corresponds to a motor neuron. (C) Twitch force of recruited fibres as a function of the total force produced in BISH muscle simulation. There is a good correspondence between this result and the data reported in [11,12]. (D) The total muscle force as a function of the percentage of motor units that have been recruited (in agreement with Walmsley et al. [13]). (E) Gradual increase in BISH muscle force during a ramped increase in input current to the monitor unit neurons. Modified from Nielsen [3].

The parameter N from eqn (1) also depends on the activity of the corresponding motor neuron. A simplified ‘synaptic-like’ model was proposed by Nielsen [3], which essentially captures the ‘catch-like’ effect [10], and was used in the present simulations to determine the frequency–response dynamics of the neuromuscular junction. Given these considerations, the total force FBISH produced by the muscle may be calculated from

 
formula
(5)

where Fsarc,i is the force produced by a typical sarcomere within the motor unit i, given the current contraction velocity of the muscle, muscle length and motor neuron activity level.

Results

The importance of modulating motor neuron firing rates is seen by plotting the activation frequency of each motor neuron as a function of the total isometric force produced by the simulated muscle (shown in Figure 1B).

The twitch force of recruited fibres as a function of the total force produced in muscles should be approximately linear [11,12] and is evidenced by plotting both quantities on a log–log plot, as shown in Figure 1(C).

A direct consequence of the muscle twitch force distribution in combination with the ordered recruitment of motor neurons is that the total muscle force will increase only slightly (to 20%) until more than half of the motor units are recruited [13,14]. This can be seen in Figure 1(D), which is a plot of the instantaneous muscle force at the moment when a given motor becomes recruited. Even when 100% of the motor units have been recruited, the BISH muscle is still only halfway towards its maximal isometric force of approx. 90 N. By continuing the ramped increase in afference to the motor neurons, the modelled muscle eventually does reach the maximal isometric tension, as shown in Figure 1(E). On the basis of similar results, Milner-Brown et al. [15] concluded that motor unit recruitment plays a role only at low levels of muscle force, whereas the increased firing rate is the main mechanism at higher force levels.

Conclusion

The importance of scaling factors in the context of muscle modelling is emphasized and a model is derived that provides quantitatively accurate results even for simulations with only a few motor units. The model accounts for many of the experimentally observable properties of motor units at the molecular and cellular levels and during whole muscle activations (see Figure 1).

Lipids, Rafts and Traffic: A Focus Topic at BioScience2004, held at SECC Glasgow, U.K., 18–22 July 2004. Edited by G. Banting (Bristol, U.K.), N. Bulleid (Manchester, U.K.), C. Connolly (Dundee, U.K.), S. High (Manchester, U.K.) and K. Okkenhaug (Babraham Institute, Cambridge, U.K.)

Abbreviations

     
  • BISH

    biceps caput breve

Part of this work was funded by the Danish National Research Foundation.

References

References
1
Huxley
A.F.
Simmons
R.M.
Nature (London)
1971
, vol. 
233
 (pg. 
533
-
538
)
2
Nielsen
B.G.
J. Theor. Biol.
2002
, vol. 
219
 (pg. 
99
-
119
)
3
Nielsen
B.G.
Ph.D. Thesis
2001
Lyngby, Denmark
Technical University of Denmark
4
Kellermayer
M.S.Z.
Smith
S.B.
Granzier
H.L.
Bustamante
C.
Science
1997
, vol. 
276
 (pg. 
1112
-
1116
)
5
Oberhauser
A.F.
Marszalek
P.E.
Erickson
H.P.
Fernandez
J.M.
Nature (London)
1998
, vol. 
393
 (pg. 
181
-
185
)
6
Nielsen
B.G.
J. Phys. Condensed Matter
2003
, vol. 
15
 (pg. 
S1759
-
S1765
)
7
Edman
K.A.P.
J. Physiol. (Cambridge, U.K.)
1988
, vol. 
404
 (pg. 
301
-
321
)
8
Hill
A.V.
Proc. Roy. Soc. B
1938
, vol. 
126
 (pg. 
136
-
195
)
9
Wilson
H.R.
J. Theor. Biol.
1999
, vol. 
200
 (pg. 
375
-
388
)
10
Burke
R.E.
Rudomin
P.
Zajac
F.E.
Brain Res.
1976
, vol. 
109
 (pg. 
515
-
529
)
11
Milner-Brown
H.S.
Stein
R.B.
Yemm
R.
J. Physiol. (Cambridge, U.K.)
1973
, vol. 
230
 (pg. 
358
-
370
)
12
Monster
A.W.
Chan
H.
J. Neurophysiol.
1977
, vol. 
40
 (pg. 
1432
-
1443
)
13
Walmsley
B.
Hodgson
J.A.
Burke
R.E.
J. Neurophysiol.
1978
, vol. 
41
 (pg. 
1203
-
1216
)
14
Purves
D.
Augustine
G.J.
Fitzpatrick
D.
Katz
L.C.
LaMantia
A.S.
McNamara
J.O.
Neuroscience
1997
Sunderland, MA
Sinauer Associates
15
Milner-Brown
H.S.
Stein
R.B.
Yemm
R.
J. Physiol. (Cambridge, U.K.)
1973
, vol. 
230
 (pg. 
371
-
390
)