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.
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  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  to find the velocity-dependent sarcomeric force Fsarc by deriving the following set of equations:
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  for details).
In , 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 , a linear combination of these may account for the double hyperbolic force–velocity relation observed by Edman .
Regardless of exactly how the sarcomeric tension is calculated (e.g. by following Hill , Huxley and Simmons  or Nielsen ), 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  for regular spiking neurons. As suggested by Nielsen , 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
where carea is an arbitrary scaling constant and ri is the radius of the ith neuron. The neuron radii were distributed according to
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
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).
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 , which essentially captures the ‘catch-like’ effect , 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
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.
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.  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.
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.)