Comparison between two muscle models under dynamic conditions

Computers in Biology and Medicine 35 (2005) 373 – 387 http://www.intl.elsevierhealth.com/journals/cobm Comparison between two muscle models under dynamic conditions R.T. Raikova∗ , H.Ts. Aladjov Bulgarian Academy of Sciences, Centre of Biomedical Engineering, Acad. G. Bonchev Str., Bl. 105, 207A, So a 1113, Bulgaria Received 4 November 2003; accepted 1 March 2004 Abstract One fundamental problem when trying to calculate the force developed by one muscle during a motor task is the muscle model. Usually, one control signal is juxtaposed to one musclotendon unit. The question is how is this signal connected to the activation of the motor units (MUs) that compose the muscle and re di erently. The aim of the paper is to compare a Hill-type muscle model to a model composed of MUs. A fast elbow exion performed by only one muscle is considered. The activation necessary for performing the motion and the corresponding frequencies are calculated for cases of fast and slow muscles using Hill-type model. Then the muscle is modelled as a mixture of 774 MUs with uniformly distributed twitch parameters. Using MotCo software the moments of impulsation of all MUs and their mechanical responses are predicted. The activation characteristics obtained by the two muscle models are compared. It is concluded that there are two essential parameters for proper muscle modelling: the lead-time and the MUs composition. ? 2004 Elsevier Ltd. All rights reserved. Keywords: Motor control; Muscle model; Motor units; Hill-model; Muscle forces; Genetic algorithm; Neural control 1. Introduction When modelling the control of the human limb motions, the nal aim is to estimate the individual muscle forces. The scienti c problem is how shall the necessary joint moments be satised by the synchronous contribution of so many muscles. A fundamental issue here is the muscle model [1]. The muscle force is usually considered a simple design variable when the external joint ∗ Corresponding author. Tel.: +3592-80-05-27; fax: +3592-72-37-87. E-mail address: [email protected] (R.T. Raikova). 0010-4825/$ - see front matter ? 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.compbiomed.2004.03.001 374 R.T. Raikova, H.Ts. Aladjov / Computers in Biology and Medicine 35 (2005) 373 – 387 Nomenclature ’ Mext d F Fmax l V FL(l) FLs FLf FV (l; V ) FVs FVf FPE1 (l) and FPE2 (l) FPE1s FPE2s FPE1f FPE2f alfa-fast alfa-slow alfa-motco alfa-motco with Tlead alfa-motco-fast f le Y f-fast f-slow fef-fast fef-slow f-motco f-motco-fast exion angle in the elbow external moment in the joint lever arm of a muscle equivalent muscle force maximal muscle force current length of the muscle normalized to the optimal length contraction velocity of the muscle force–length relationship calculated force–length relationship during the motion if the muscle is slow calculated force–length relationship during the motion if the muscle is fast force–velocity relationship calculated force–velocity relationship during the motion if the muscle is slow calculated force–velocity relationship during the motion if the muscle is fast forces in the passive elements calculated forces in the passive elements if the muscle is slow calculated forces in the passive elements if the muscle is fast muscle activation calculated activation in case of fast muscle calculated activation in case of slow muscle “activation parameter” calculated by MotCo software in case of uniformly distributed parameters of the motor units “activation parameter” calculated by MotCo software in case of uniformly distributed parameters of the motor units but shifted forward with 100 ms “activation parameter” calculated by MotCo software in case when all motor units are fast-like stimulus frequency muscle e ective length yielding calculated frequency during the motion in case of fast muscle calculated frequency during the motion in case of slow muscle calculated frequency during the motion in case of fast muscle accounting for the yielding and using the muscle e ective length calculated frequency during the motion in case of slow muscle accounting for the yielding and using the muscle e ective length average frequency of the motor units ring calculated by MotCo software in case of uniformly distributed parameters of the motor units average frequency of the motor units ring calculated by MotCo software in the case when all motor units are fast-like R.T. Raikova, H.Ts. Aladjov / Computers in Biology and Medicine 35 (2005) 373 – 387 375 moments, calculated by using inverse dynamic approach, are distributed among the muscles driving the correspondent joints on the base of some optimization criterion [2]. To validate the predicted muscle forces, usually processed electromyographic signals (EMGs) are used. That is why, the question about the relationship between EMGs and correspondent muscle forces arises naturally [3]. Other approaches are based on Hill-type muscle models. The constitutive equations describing this model behaviour require data for current length and contraction velocity of the muscle. Knowing the latest, the individual muscle forces can be computed if the so-called “activation” parameters [4] are known as functions of time. They can be calculated, again using optimization, but if the design variables are the activation parameters [5–8]. Another possibility is to use directly processed EMGs, i.e. to estimate the muscle activation on the basis of EMG signals [9] and using Hill-type model to calculate the muscle forces. Then summing the moments of all muscles, the total joint moment can be obtained and compared to an experimentally measured joint moment [10]. A critical point is the EMGs that is necessarily used as an estimation of the muscle activation or “neural input” [11]. A common paradigm for all these approaches is “one muscle— one control signal” [12]. The muscle force, however, is graduated by including/excluding motor units (MUs) or by changing their discharge pattern, i.e. through recruitment and rate coding [13]. Thus, the muscle force is controlled not by one smooth, continuous signal, but through discrete impulses from motoneurons. How can they be related to the surface EMGs and to the so-called “activation” parameter in Hill-type muscle models, is a question of great importance. The ring of all MUs in a muscle during a motion is impossible to record in vivo. The same refers to the forces of the individual MUs. Modelling tools can help much in respect to MUs action. The recently proposed model in [14] considers the muscle as a mixture of many MUs with di erent peculiarities. A hierarchical genetic algorithm that mimics the action of the nervous system is implemented. It chooses the time moments of impulsation of all MUs so that a preliminary given motion to be performed. Since the task is highly undetermined, di erent time-dependent criteria can be imposed simultaneously on the solution. Such criteria can be: the minimal deviation of the calculated joint moment from the desired one, the minimal muscle force, the minimal joint reactions, the minimal activation and so on. The developed MotCo software (http://motco.dir.bg) allows us to obtain di erent characteristics of the modelled muscles. These characteristics can be of mechanical nature (like total muscle force, forces of all MUs during the motion) and of activation nature (like time moments of impulsation of all MUs, interpulse intervals, number of impulses per given time period) [15]). That is why it seems to be a suitable tool for detailed investigation of the relationship between the discrete impulsation of the MUs and a continuous signal re ecting the whole muscle activity. The previously reported simulations [5,14] were directed to investigations of the force-sharing problem and it was supposed that three exors and two extensors perform the simulated elbow motions. In this case many parameters in uence the predicted results, and this especially refers to the optimization criterion for muscle forces’ distribution. Considering one degree-of-freedom model with a single muscle is useful when the aim is (1) to clarify the nature of the parameter “activation” used in Hill-type muscle models and (2) to study its relationship with MUs functioning. The aim of the paper is to compare a Hill-type muscle model and the model based on MUs simulating a fast elbow exion in the sagittal plane performed by one muscle only. The study is directed to activation characteristics and dependence of the simulation results on the muscle composition of MUs with di erent peculiarities. 376 R.T. Raikova, H.Ts. Aladjov / Computers in Biology and Medicine 35 (2005) 373 – 387 2. Simulated motion and the two muscle models A fast elbow exion (from fully extended forearm, ’ = 0◦ , to fully exed forearm, ’ = 150◦ ) with duration of 0:35 s in the sagittal plane has been simulated. It is performed by one muscle only, namely biceps brachii (BIC), as a muscle equivalent. The external joint moment (Mext ) during the motion is calculated in accordance with a sinusoidal law for angle displacement (for details, see [14]). The muscle lever arm (d) depends on the joint angle (’) and is calculated using literature regression equations [16]. 2.1. Hill-type model The used Hill-type model is strictly based on the equations given in Brown et al. [17] and the constants given in Table 1, p. 640 in the same paper [17]. The main equation that relates the muscle force F to the unknown activation (it is de ned as an “activation–frequency relationship” in [17] and represents a measure of “the fraction of available cross-bridges that are cycling in all recruited muscle bres”) is F (1) = [FL(l)FV (l; V ) + FPE2 (l)] + FPE1 (l); Fmax where 0 6 6 1, FPE1 (l) and FPE2 (l) represent passive elastic elements, FL(l) is the force–length relationship, FV (l; V ) is the force–velocity relationship, l and V are the current length and contraction velocity of the muscle (normalized to the optimal length at which the muscle develops maximal isometric force), Fmax is the maximal possible force of the muscle, being calculated by multiplying the physiological cross-sectional area of the muscle, namely 5:37 cm2 , by 100 N=cm2 . 1 The aim is to nd the muscle activation and the corresponding frequency that can generate it, so that the motion to be performed. Since the muscle force necessary for motion performing has to satisfy the moment equation Mext = Fd, can be calculated from Eq. (1) as (Mext =dFmax ) − FPE1 (l) (2) : = FL(l)FV (l; V ) + FPE2 (l) According to [17] the activation–frequency relationship, i.e. (f), can be described by a function like 
1 −(f=a1 n)n (3) −1 ; (f) = 1 − e ; n = a2 + a 3 ‘ where a1 , a2 and a3 are constants. Developing further this equation according to experimental data, Brown et al. [17] introduces terms “e ective length”, le , (time-delay version of length accounting for the previous state of the muscle, it depends on l) and “yielding”, Y (that depends on V ). Then the activation–frequency relationship becomes (for more details see in [17]) 
1 −(Yf=b1 p)p (4) −1 ; (f) = 1 − e ; p = b2 + b 3 ‘e where b1 , b2 and b3 are constants that have di erent values for fast and slow MUs. 1 The upper limit of this value (see in [1]) is used here since the modelled muscle represents a muscle equivalent and involves the action of other elbow exors as brachialis and brachioradialis. R.T. Raikova, H.Ts. Aladjov / Computers in Biology and Medicine 35 (2005) 373 – 387 377 Calculating from Eq. (2), Eqs. (3) and (4) can be solved with respect to the frequency f. Note that this frequency is also normalized to this one that is necessary to produce one-half of the maximal muscle force at the optimal length. In Brown et al. [17] all constants are given separately for fast- and slow-twitch-muscles, only FPE2 is one and the same for both cases. Human limb muscles, however, are composed of di erent portions of slow and fast MUs [18]. Since the activation parameter cannot be divided into fast and slow parts, without using optimization tools, then two distinct cases are considered here: (1) the muscle consists of fast MUs only, i.e. it is a fast muscle and (2) the muscle consists of slow MUs only, i.e. it is a slow muscle. 2.2. Muscle model based on MUs Three simulations are performed using MotCo software. (1) The muscle consists of 774 MUs, that are not separated into slow or fast MUs, but their parameters are uniformly distributed from the fastest MU to the slowest MU. It means that the distribution of the maximal amplitude of the MUs twitches is from 0.5587 to 0.06491 N, the contraction time varies from 30 to 80 ms, the half-relaxation time changes from 60 to 200 ms, and the lead time varies from 20 to 70 ms. (2) The muscle consists of 774 fast-like MUs and their parameters are again uniformly distributed, but within other limits: the maximal amplitude of the individual twitches is from 0.399 to 0:2318 N, the contraction time changes from 30 to 40 ms, the half-relaxation time varies from 60 to 100 ms, and the lead time varies from 20 to 30 ms. (3) The muscle consists of 774 slow-like MUs and their parameters are again uniformly distributed, but within other limits: the maximal amplitude of the individual twitch is from 0.3601 to 0:1801 N, the contraction time is from 70 to 80 ms, the half-relaxation time is from 175 to 200 ms, the lead time is from 60 to 70 ms. The developed hierarchical genetic algorithm simulates the action of the nervous system and chooses the time moments of impulsation of all MUs, so that the motion to be performed. Knowing the impulsation of all MUs and summing their mechanical responses, the total muscle force F and the correspondent joint moment are calculated as functions of time. The optimization criteria currently included in the tness function [5], which estimate the current solution’s reliability, are: minimal deviation of the calculated joint moment (obtained by predicted MUs forces) from the required joint moment Mext with a weight coecient of 100, minimal total muscle force with a weight of 1, and minimal number of impulses of all MUs with a weight of 1 during the whole motion. 3. Results Parameters of the simulated motion (joint angle, velocity and acceleration) as functions of time are presented in Fig. 1. The muscle BIC lever arm, d, the external joint moment Mext and the muscle force F necessary for satisfying this moment (that is calculated as F = Mext =d) are also shown. The current length l and the contraction velocity V of the muscle during the motion (normalized 378 R.T. Raikova, H.Ts. Aladjov / Computers in Biology and Medicine 35 (2005) 373 – 387 16 12 8 4 0 0 -4 -8 -12 -16 0.05 0.1 angle velocity acceleration/10 Mext d F/10 0.15 0.2 0.25 0.3 0.35 t [s] Fig. 1. Parameters of the modelled motion: joint angle in radians, angle velocity in radians per second, angle acceleration in radians per second 2 , the necessary joint moment Mext in N m, the muscle lever arm, d, in cm, the calculated muscle force, F, in N. For a better view the true acceleration and muscle force are divided to a factor of 10. to the optimal muscle length) are calculated from the regression equations given in [16,19] (see Fig. 2). Using the equations given in Brown et al. [17] the forces of the passive and contractile elements are calculated separately for slow and fast muscles. FPE2 was equal to zero during the whole motion. FPE1 was near zero for both slow and fast muscles (i.e. FPE1s and FPE1f in Fig. 2). Hence, the in uence of the passive elastic elements can probably be neglected for such dynamic conditions, and the muscle force can be approximated as F ≈ FL(t)FV (t) [20,21]. The curve FL(t) does not di er signi cantly for slow (FLs ) and fast (FLf ) muscles (see Fig. 2), but the curve FV (t) for slow muscle (FVs in Fig. 2) is enough under the one for fast muscle (FVf in Fig. 2). Hence, the capabilities of the slow MUs for force developing are much less when the contraction velocity is high. The activation parameters calculated by using Eq. (2), namely alfa-fast in case of fast muscle and alfa-slow in case of slow muscle (see Fig. 3), have a similar form to the form of the curve F=Fmax , but the amplitude of alfa-slow is much bigger. A similar parameter, alfa-motco is derived from the rst simulation with MotCo software. It is calculated from the time moments of impulsation of all MUs (Fig. 4a). In order to show how is the total muscle force formed from the action of the separate MUs, the mechanical responses of the correspondent MUs to this activation are shown in Fig. 4b. For better view only the rst 194 (the fastest) MUs are shown. Summing the forces of all individual MUs (Fig. 4b), the required muscle force F (shown in Fig. 1) is reached and the necessary joint moment Mext is satis ed with great accuracy. Alfa-motco is calculated as follows: for a xed time interval (the currently used interval is 120 ms, it is equal to 2.2 times the average contraction time for all muscle MUs and has been experimentally selected) the number of received impulses from one MU is counted and divided to the maximal possible number of impulses to R.T. Raikova, H.Ts. Aladjov / Computers in Biology and Medicine 35 (2005) 373 – 387 379 1.5 1 0.5 0 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 t [s] -0.5 -1 -1.5 -2 l V FPE1f FPE1s FLf FLs FVf FVs Fig. 2. Parameters for the Hill-type muscle model: length l and contraction velocity V of the muscle during the motion, force of the passive element in case of fast (FPE1f ) and slow (FPE1s ) muscle, force–length (FLf ; FLs ) and force–velocity (FVf ; FVs ) relationships in case of fast and slow muscle, respectively. All curves are normalized to the optimal muscle length that is chosen to be 36 cm. The maximal length of the muscle is 37:806 cm, the minimal is 28:466 cm [9,19]. which this MU can respond (this number depends on the absolute refractory period of the current MU, which is supposed equal to the contraction time). This time interval is applied as a moving frame, the obtained parameter being averaged for all MUs. Since the activation calculated by using Eq. (2) does not take into consideration the electro mechanical delay [22], with a purpose of better comparison the parameter alfa-motco is shifted forward with a time period of 100 ms (i.e. alfa-motco with Tlead). As it can be seen from Fig. 3, this shifted curve has similar form to the one obtained by Eq. (2) but its amplitude is much smaller. The activation parameter calculated for the case when all muscle MUs are fast-like, i.e. for the second simulation, is also shown in Fig. 3 (the line noted as alfa-motco-fast). When all muscle MUs were slow-like, i.e. the third simulation, it was unable to satisfy the joint moment (see discussion) and that is why the correspondent activation parameter misses in the gure. The frequencies calculated by using Eq. (3) (i.e. f-slow in case of slow muscle and f-fast in case of fast muscle) and Eq. (4) (i.e. fef-slow in case of slow muscle and fef-fast in case of fast muscle) are shown in Fig. 5. They are multiplied to a value of 34 Hz, i.e. to the normalized frequency used in Brown et al. [17]. There are no signi cant di erences between these two formulas in case of fast muscle (i.e. f-fast versus fef-fast), but if the muscle is slow, the frequencies calculated using the parameters “yielding” and “e ective length” (i.e. Eq. (4)) and using the simple Eq. (3), are 380 R.T. Raikova, H.Ts. Aladjov / Computers in Biology and Medicine 35 (2005) 373 – 387 0.8 F/Fmax alfa-fast alfa-slow alfa-motco alfa-motco with Tlead alfa-motco-fast 0.6 0.4 0.2 0 -0.15 -0.05 0.05 0.15 0.25 0.35 0.45 t [s] Fig. 3. Calculated activation versus the normalized muscle force F=Fmax . alfa-fast—activation calculated using Eq. (2) in case of fast muscle; alfa-slow–activation calculated using Eq. (2) in case of slow muscle; alfa-motco—activation parameter calculated by using MotCo software in case when the MUs parameters are uniformly distributed; alfa-motco with Tlead —the same curve as alfa-motco but shifted forward with 100 ms; alfa-motco-fast—activation parameter calculated by using MotCo software in case when the muscle MUs are fast-like. essential di erent (i.e. fef-slow versus f-slow). As it can be seen, a slow muscle has to be stimulated with fairly bigger frequency than a fast one in order to have the same activation parameter and to develop the same force (i.e. f-slow versus f-fast or fef-slow versus fef-fast). From the predicted time moments of impulsation of the MUs with MotCo software (Fig. 4a), we calculated an average frequency of ring of MUs of our modelled muscle. It is presented with the line f-motco in Fig. 5 for the case when the muscle has mixed MUs (i.e. the rst simulation). This parameter presents the number of impulses received by all MUs for a xed time period of 1 ms (the so-called “simulated EMG” [5]). In order to obtain an average frequency in Hertz this number is divided to the number of MUs, namely 774, and multiplied by 1000. The temporal frequency curve for the case when all MUs were fast-like, i.e. f-motco-fast in Fig. 5, was similar to the curve f-motco, but its amplitude was slightly higher. Up to now only the results from the rst two simulations with MotCo software have been shown. Considering the third case, when all MUs were slow-like the software was unable to satisfy the joint moment by normal impulsation of the MUs (Fig. 6b) opposite to the rst two cases (Fig. 6a). The satisfaction of the joint moment in the beginning of the motion leads to the impossibility to follow its fast decrease phase (Fig. 6c)—see discussion. R.T. Raikova, H.Ts. Aladjov / Computers in Biology and Medicine 35 (2005) 373 – 387 381 Fig. 4. Time moments of impulsation and mechanical responses of the rst 194 MUs (the most fast) predicted by MotCo software. Fig. 4a. Each point corresponds to one neural impulse, each row represents one MU, the bars in the right correspond to the total number of impulses received by the respective MU (the maximal value is 5). Note that the time scale is larger than the duration of the motion because of the lead time and the total duration of a twitch. Fig. 4b. The mechanical responses of the rst 194 MUs to the correspondent activation given in Fig. 4a. 382 R.T. Raikova, H.Ts. Aladjov / Computers in Biology and Medicine 35 (2005) 373 – 387 f-slow 60 fef-slow f-fast fef-fast 50 f-motco f [Hz] f-motco-fast 40 30 20 10 f-motco 0 -0.1 0 0.1 0.2 0.3 0.4 t [s] f-motco-fast 0 -0.1 0 0.1 0.2 0.3 0.4 t [s] Fig. 5. The computed frequencies. f-slow and f-fast are calculated using Eq. (3) and considering that the muscle is slow or fast, respectively; fef-slow and fef-fast are calculated using Eq. (4) considering that the muscle is slow or fast, respectively (all values are multiplied by 34 Hz, i.e. the normalized frequency used in Braun et al. [17]); f-motco—an average frequency calculated by MotCo software for the case of uniformly distributed MUs (it presents the number of impulses applied on all MUs of the muscle for each one millisecond time interval, this number is then divided to the total number of the MUs, namely 774); f-motco-fast—an average frequency calculated by MotCo software for case of fast-like MUs. 4. Discussion For many investigation and application purposes the possibility to derive one control signal for a muscle during its working proved to be very attractive. It is questionable whether this is possible to be well done because of the discrete nature of the impulsation of the MUs, i.e. of their ring. Moreover, a new hypothesis was presented in [23] that the muscle force is not directly related to the activation. Since the EMGs is main accessible evidence of the muscle state, its relationship with this control signal is of great importance. This problem was considered in the present paper for dynamic conditions. It was supposed that a fast elbow exion is performed by one muscle only, aiming to avoid the in uence of the optimization criteria when there is a synchronous action of many synergistic and antagonistic muscles. R.T. Raikova, H.Ts. Aladjov / Computers in Biology and Medicine 35 (2005) 373 – 387 383 Fig. 6. Desired (bold line) and calculated (grey line) joint moment using MotCo software. (a) the muscle is composed of fast-like MUs (see Section 2.2 for MUs parameters); (b) the muscle is composed of slow-like MUs and the required joint moment cannot be well described by summing the forces of these MUs through normal impulsation; (c) the muscle is composed of slow-like MUs but the conditions during the simulation are changed—the genetic algorithm is forced to satisfy only the acceleration phase of the movement, i.e. from 0 to 120 ms. Di erent Hill-type muscle models are reported in the literature, but a few authors mention exactly which constants are used in the constitutive equations and what are their assumptions about the percent of fast and slow muscle parts. That is why we used the equations reported in Brown et al. [17], where these constants are obtained by natural animal experiments separately for fast and slow muscles. Our results showed that for dynamic conditions the activation parameter depended rather on the type of the muscle (Fig. 3)—the maximal value of alfa-fast is about 2.7 times lower than the maximal value of alfa-slow. The main reason is the di erence in the curves “force–velocity” for fast and slow muscles (Fig. 2—FVf versus FVs ), while the “force–length” relationship is no much sensitive to the muscle type (FLf versus FLs ). The activation parameter calculated by using the ring of the MUs predicted by MotCo software (Fig. 4a) was even smaller than alfa-fast (see Fig. 3). Shifting it forward with 100 ms, it can be seen that all activation parameters have forms similar to the curve F=Fmax , but di erent amplitudes. Surprisingly, for the case when all MUs were faster, 384 R.T. Raikova, H.Ts. Aladjov / Computers in Biology and Medicine 35 (2005) 373 – 387 the activation increased in amplitude (i.e. alfa-motco-fast versus alfa-motco) contrary to the results obtained by Hill-type muscle (i.e. alfa-fast versus alfa-slow). This probably is due to the chosen lower maximal force of the twitches of the fast-like MUs, namely 0.399 versus 0:5587 N for the rst case with uniformly distributed MUs. That is why for reaching the same muscle force with weaker MUs they have to be red with higher frequency. We shall also emphasize that the time delay between alfa-motco-fast and the muscle force is smaller in comparison with alfa-motco. The curve alfa-motco-fast has to be shifted forward with 70 ms in order to account for electro-mechanical delay. Therefore, the activation parameter depends strongly on the muscle’s composition of di erent MUs. It can change if the time parameters of the MUs twitches (contraction time, lead time, half-relaxation time) are changed, as well if the maximal twitch forces are di erent. The time delay between activation and force development is not accounted in Eq. (2). It can be accounted during the process of converting the activation in EMGs (or vice versa), i.e. in the excitation–contraction dynamics. The presented simulations with MotCo software showed that this time delay depends on the type of the MUs composing the muscles and it was between 100 and 70 ms (according to Hof [3] it varies from 50 to 200 ms). Having considered the muscle a composition of slow-like MUs only, the genetic algorithm could not generate suitable impulsation of the MUs so that the required joint moment to be satis ed (Fig. 6b) since the motion was quite fast. The values used for contraction and half-relaxation times of these MUs are quite large, but they are in accordance with the data reported in the literature [24]. If the initial impulsation is arti cially forced so that the fast initial increase of the joint moment be satis ed, then the residual MUs mechanical activity decreases much slowly than the joint moment (Fig. 6c). If an antagonistic muscle was presented in the model its action could help to perform the motion since its moment would be with negative sign and its action could compensate for the residual positive moment after the peak. The Hill-type models have no such restrictions, their only limit is 0 6 6 1. That is why the results they predict can be unrealistic from the physiological point of view. For example, the maximal frequency calculated by using Eq. (4) is 62:4 Hz, which is probably impossible for slow human MUs, since it presumes an interpulse interval of about 15 ms. The slow-like MUs in our model can be red with maximal frequency of 14:3 Hz. Of course, the maximal values of the frequencies shown in Fig. 5 and calculated by Eq. (3) and (4), namely f-slow, f-fast, fef-slow and fef-fast, depend on the normalizing value of 34 Hz, reported in [17], but obviously some physiological limits must be imposed using Hill-type muscle models. The average frequencies calculated by MotCo software, f-motco and f-motco-fast, are not such raw, smooth curves as the correspondent curves for Hill-type model (Fig. 5). They look like an EMG signal and reach the maximal value of 16:79 Hz while for fef-slow this maximum is 62:4 Hz. Our opinion is that this output of MotCo software is more realistic since the ring frequency of MUs is not constant during the motion, while the Eqs. (3) and (4) are obtained using experimental, arti cial but not natural, stimulation with constant frequencies (see also the discussion in [25] where these equations are elaborated). There are other approaches to relate the activation parameter with EMGs (or vice versa) [7,12,4,26,27], i.e. activation–excitation relationship. They were not used in the present paper because of the inability to nd a paper where the necessary constants for this relationship are given separately for fast and slow muscles. Some calculations were performed, however, using the equation given in [4, p. 371], that relates muscle activation with neural excitation. No matter how many di erent constants have been tried, the curve of the excitation was very close to this of the activation and it was not time shifted, i.e. the electro-mechanical delay was not accounted again. R.T. Raikova, H.Ts. Aladjov / Computers in Biology and Medicine 35 (2005) 373 – 387 385 Both muscle models (Hill-type and MUs based) used in this paper have their sources of error [21]. The most critical factors in the model based on MUs are the parameters of the twitches. The rst simulated case, where these parameters were uniformly distributed, seems more realistic since there are evidences that the contractile properties of human limbs’ MUs span a wide range [24], but are not strictly separated to fast and slow ones. 5. Summary The purpose of this paper was to compare a Hill-type muscle model based on the work of Brown et al. [17] and the new muscle model based on motor units [14], which have been developed by the authors. Attention was paid to the net activation parameter that re ects the whole muscle activity and its relationship with the MUs ring. Supposing that a fast elbow exion is performed by only one muscle, and calculating the muscle force necessary for satisfying the joint moment, the activation and the frequencies were obtained using the Hill-type model. They were di erent for cases of fast and slow muscle and did not account for the electro-mechanical delay. It was also observed that due to a lack of suitable constraints such model can predict results that are impossible in a real situation (for example performing a very fast motion from a muscle composed of very slow MUs). Hence, a big attention had to be paid (1) on the constants used in the constitutive equations for describing the relationships force–length–velocity–activation (they are quite di erent for fast and slow muscles) and (2) on the way the surface EMGs is related to the predicted muscle force, activation or frequency, especially on the electro-mechanical delay. The model developed by the authors and the MotCo software allow the muscle’s simulation as if composed of a real number of di erent types of MUs. Thus, the relationship between the action of individual MUs (their ring and mechanical activity) and parameters related to the whole muscle can be studied in details. The net activation parameter that is an output of the software, had similar form to the activation curve obtained by the Hill-type model, but the electro-mechanical delay was accounted properly. The other output, the so-called “simulated EMG” (de ned as an average frequency of the MUs ring) is much closer to an EMG signal, than the frequencies obtained using the equations given in Brown et al. [7]. In our opinion it is more realistic since the MUs do not re with constant and continuous frequencies during a motion. These software outputs depend on the type of MUs composing the muscle. If there are not enough suitable MUs, the required motion cannot be performed. A further comparison between experimental and simulation results will help a lot in the practical application of this new model and software. Appendix. c c FL(l) = e−|(l 1 −1)=c2 | 3  c4 − V    c + (c + c l)V 4 5 6 FV (l; V ) = 2    c7 − (c8 + c9 l + c10 l )V c7 + V if V 6 0; if V ¿ 0; 386 R.T. Raikova, H.Ts. Aladjov / Computers in Biology and Medicine 35 (2005) 373 – 387 FPE1 (l) = c11 ln(e(l−c12 )=c13 + 1) FPE2 (l) = c14 [ec15 (l−c16 ) − 1] if FPE2 (l) ¡ 0 else FPE2 (l) = 0 ci (i = 1; 2; : : : ; 14) are constants (see Table 1, [17, p. 640]). References [1] A.J. van den Bogert, K.G.M. Gerritsen, G.K. Cole, Human muscle modelling from a user’s perspective, J. Electromyog. Kinesiol. 8 (1998) 119–124. [2] D. Tsirakos, V. Baltzopoulos, R. Bartlett, Inverse optimization: functional and physiological considerations related to the force-sharing problem, Crit. Rev. Biomed. Eng. 25 (1997) 371–407. [3] A.L. Hof, The relationship between electromyogram and muscle force, Sportverletzung Sportschaden 11 (1997) 79–86. [4] F.E. Zajac, Muscle and tendon: properties, models, scaling, and application to biomechanics and motor control, Crit. Rev. Biomed. Eng. 17 (1989) 359–410. [5] R. Raikova, H. Aladjov, The in uence of the way the muscle force is modeled on the predicted results obtained by solving indeterminate problems for a fast elbow exion, Comput. Methods Biomech. Biomed. Eng. 6 (2003) 181–196. [6] H. Rehbinder, C. Martin, A control theoretic model of the forearm, J. Biomech. 34 (2001) 741–748. [7] D.J. Thelen, F.C. Anderson, S.L. Delp, Generating dynamic simulations of movement using computed muscle control, J. Biomech. 36 (2003) 321–328. [8] G. Li, K. Kawamura, P. Barrance, E.Y.S. Chao, K. Kaufman, Prediction of muscle recruitment and its e ects on joint reaction forces during knee exercises, Ann. Biomed. Eng. 26 (1998) 725–733. [9] L. Wang, T.S. Buchanan, Prediction of joint moments using a neural network model of muscle activations from EMG signals, IEEE Trans. Neural Systems Rehabilitation Eng. 10 (2002) 30–37. [10] K. Manal, R.V. Gonzalez, D.G. Lloyd, T.S. Buchanan, A real-time EMG-driven virtual arm, Comput. Biol. Med. 32 (2002) 25–26. [11] J.P. van Zandwijk, M.F. Bobbert, J. Harlaar, Predictions of mechanical output of the human m.triceps surae on the basis of electromyographic signals: the role of stimulation dynamics, J. Biomech. Eng. 122 (2000) 380–386. [12] K. Manal, T.S. Buchanan, A one-parameter neural activation to muscle activation model: estimating isometric joint moments from electromyograms, J. Biomech. 36 (2003) 1197–1202. [13] R.B. Stein, Peripheral control of movement, Physiol. Rev. 54 (1974) 215–243. [14] R.T. Raikova, H.Ts. Aladjov, Hierarchical genetic algorithm versus static optimization—investigation of elbow exion and extension movements, J. Biomech. 35 (2002) 1123–1135. [15] R.T. Raikova, Hr.Ts. Aladjov, Simulation of the motor units control during a fast elbow exion in the sagittal plane, J. Electromyog. Kinesiol. 14 (2004) 227–238. [16] M.A. Lemay, P.E. Crago, A dynamic model for simulating movements of the elbow, forearm, and wrist, J. Biomech. 29 (1996) 1319–1330. [17] I.E. Brown, E.J. Cheng, G.E. Loeb, Measured and modelled properties of mammalian skeletal muscle. II. The e ects of stimulus frequency on force–length and force–velocity relationships, J. Muscle Res. Cell Motility 20 (1999) 627–643. [18] M.A. Johnson, J. Polgar, D. Weightman, D. Appleton, Data on the distribution of ber types in thirty six human muscles. An autopsy study, J. Neurol. Sci. 18 (1973) 111–129. [19] P. Pigeon, L.H. Yahia, A.G. Feldman, Moment arms and lengths of human upper limb muscles as function of joint angles, J. Biomech. 29 (1996) 1365–1370. [20] B.M. van Bolhuis, C.C.A.M. Gielen, A comparison of models explaining muscle activation patterns for isometric contractions, Biol. Cybern. 81 (1999) 249–261. [21] E.J. Perreault, C.J. Heckman, T.G. Sandercock, Hill muscle model errors during movements are greatest within the physiologically relevant range of motor unit ring rates, J. Biomech. 36 (2003) 211–218. R.T. Raikova, H.Ts. Aladjov / Computers in Biology and Medicine 35 (2005) 373 – 387 387 [22] D.A. Gabriel, J.P. Boucher, E ects of repetitive dynamic contractions upon electromechanical delay, J. Appl. Physiol. 79 (1998) 37–40. [23] A.L. Hof, Muscle mechanics and neuromuscular control, J. Biomech. 36 (2003) 1031–1038. [24] K.M. Chan, T.J. Doherty, W.F. Brown, Contractile properties of human motor units in health, ageing, and disease, Muscle Nerve 24 (2000) 1113–1133. [25] I.E. Brown, G.E. Loeb, Measured and modeled properties of mammalian skeletal muscle: IV. Dynamics of activation and deactivation, J. Muscle Res. Cell Motility 21 (2000) 33–47. [26] J.-Yu. Cheng, T.S. Buchanan, Comparison of three EMG-driven muscle models: performance for elbow joint under time-varying loading, Proceedings of ASME Summer Bioengineering Conference, 1999, pp. 559–560. [27] M.G. Pandy, F.E. Zajac, E. Sim, W.S. Levine, An optimal control model for maximum-height human jumping, J. Biomech. 23 (1990) 1185–1198. Rositsa T. Raikova received her Ph.D. degree in Biomechanics from the Institute of Mechanics and Biomechanics at the Bulgarian Academy of Sciences in 1993 and her M.S. degree from the So a University, Faculty of Mathematics and Mechanics, in 1979. She graduated in Robotics in Technical University of So a from 1980 till 1982. Now she is Associate Professor at the Centre of Biomedical Engineering, Bulgarian Academy of Sciences. Her research interests are in the eld of biomechanics and motor control of the human upper limbs. Hristo Ts. Aladjov graduated in Computer Engineering in 1998, later he received his Ph.D. degree in Arti cial Intelligence in 2002. Since 1998 he works in the Centre of Biomedical Engineering, Bulgarian Academy of Sciences as Research Assistant. His major research interests are related to the application of the methods of arti cial intelligence for modelling in biology and particularly in motor control and biomechanics. He is also interested in new algorithm development and their implementation.