Simplified Cost Functions Meet Advanced Muscle Models to Streamline Muscle Force Estimation

1. Introduction

Human motion is powered by the activation of the skeletal muscles, which respond to electrochemical signals from the nervous system. The adult human musculoskeletal system, consisting of approximately 350 joints and more than 650 named muscles, presents a complex system for sharing muscle forces [1,2]. In particular, many muscle groups control the same degree of freedom (DoF) across multiple joints, leading to a overactuated systems, where infinitely many possible muscle force distributions can generate the same net torques at the joints. Therefore, estimating muscle forces solely based on torque data at the joints poses a conceptual modeling challenge [3,4].

The brain plays a pivotal role in the activation of muscles, ensuring that they work harmoniously to produce the desired movement. A fundamental question in musculoskeletal modeling is how the brain activates these muscles to manage such an over-actuated system. Typically, researchers mimic this process through optimization strategies that involve designing cost functions. These functions can directly represent physiological processes based on motor control theories or use abstract mathematical formulas that are not directly related to physiological functions [5]. These optimization strategies help address redundancy in muscle activation, approximating a strategy to distribute the torque contributions of each muscle based on predefined assumptions [6,7].

The actual mechanisms by which the central nervous system controls muscle activation are still largely unknown. Consequently, various theoretical cost functions have been proposed to approximate these neuromuscular control strategies [8,9,10,11]. Some of these functions aim to simulate the complex biochemical pathways influencing muscle contractions [12], while others focus on minimizing mechanical energy across movements [13]. Although trends have evolved towards developing sophisticated models, the relationship between the complexity of the cost functions and the accuracy of the results is not definitive. In fact, simpler models have occasionally produced results that align closely with the experimental data [14]. Furthermore, in contexts where speed, ease of implementation, and clarity are crucial, such as educational uses, performance evaluations, or preliminary analyses, simpler, more intuitive cost functions may become necessary or at least be more advantageous.

Recent work has proposed cost functions that leverage both kinematic and electromyography (EMG) data to calculate the force distribution across muscle movements [15]. This approach effectively minimizes mechanical stress and accounts for muscle co-contractions by aligning model predictions closely with empirical EMG results. Although this methodology has advanced our understanding of muscle dynamics, it does not consider alternative cost function strategies that focus on minimizing individual muscle forces, muscle force per physiological cross-sectional area (PCSA), or mechanical work produced by each muscle [16,17,18,19]. Furthermore, the use of EMG data introduces complexities and potential inaccuracies into the model that may limit the practicality of such methods in applications where accurate muscle force estimation is less critical. Our study builds on this foundational work by examining simpler cost functions in conjunction with detailed muscle models, with the aim of enhancing both the accessibility and applicability of muscle force estimation across various domains.

In this study, we use a simple one DoF bicep curl model to explore how different muscle models and cost functions influence the distribution of muscle force. Our goal is to demonstrate that even less complex cost functions can yield highly accurate muscle force estimates when combined with detailed muscle models. We evaluated three muscle models of increasing complexity: (a) a basic force model with two contact points[6,14], (b) a Hill-type muscle model incorporating active force-length and force-velocity relationships [20], and (c) a Hill-type model that also includes passive force-length relationships [13]. We pair these models with the minimization of three different cost functions that do not require EMG data: (a) total muscle force [14,16,19], (b) total work performed by the muscles [17], and (c) total stress developed in the muscles, where stress is the force per physiological cross-sectional area (PCSA) of the muscle [9,14,18,19]. The choice of cost function is based on its suitability for the specific muscles and tasks involved, as noted in the literature [6]. For example, stress minimization is often used in arm muscle modeling [6,15,21], whereas energy expenditure minimization is preferred in leg muscles for gait analysis [6,12,22]. We evaluate the efficacy of each model and the cost function combination by comparing the computed muscle and arm moments with the established elbow flexion data. This comparison assesses the accuracy of our models and helps identify the most effective strategies to simulate realistic muscle actions.

2. Methodology

To explain our methodology, we begin by describing the biomechanical model for the bicep curl motion with one degree of freedom (DoF) (see Section 2.1), which we use to investigate the impact of different combinations of muscle activation models and optimization strategies on muscle force estimation. Next, in Section 2.2, we detail the process for calculating the forces and torque required at the elbow joint for the curling motion. We then examine various muscle models (detailed in Section 2.3) and cost functions (outlined in Section 2.4), which we employ to estimate the forces generated by the elbow flexors, taking into account an applied load and the prescribed angular velocity of the curl. Finally, Section 2.5 describes how we validated our model results.

2.1. Biomechanical Model

The biomechanical model of the arm was developed using MATLAB Simscape Multibody (version 4.16, The Mathworks, Natick, MA), as a multibody system with two linked segments representing the upper arm (humerus) and the forearm. The elbow joint, comprising the radiohumeral, ulnohumeral, and proximal radioulnar joints, primarily allows two degrees of freedom (DoF): flexion / extension (bending and unbending of the elbow) and pronation / supination (twisting the forearm about its longitudinal axis) [23]. However, the bicep curl is dominated by flexion / extension and therefore in this study the joint was approximated to have a single DoF, similar to a single-axis hinge.

The geometric and inertial parameters of the forearm were obtained from existing literature [12,14] and are summarized in Table 1 . The model includes the three primary flexor muscles of the elbow: the biceps brachii (bic), brachialis (bra) and brachioradialis (brd). The origins, insertions, and orientations of the muscles relative to the linked segments are shown in Figure 1 and detailed in Table 1 [24,25,26]. The orientation of each muscle is characterized by the parameter ${\alpha }_i$, where i refers to the muscles: bic, bra, or brd. This parameter ${\alpha }_i$ represents the acute angle between the line of action of the muscle force and the longitudinal axis of the forearm. Angle ${\alpha }_i$ determines the moment arm for each muscle, influencing the calculation and behavior of the muscle forces during motion.

In our simulation case study, we consider a range of forearm rotations, in which the angle $\theta$ varies from ${15}^{\circ }$ to ${120}^{\circ }$ (see Figure 1). The forearm is considered to rotate at a constant angular velocity $\omega$ while carrying a mass W = $30\mathrm{N}$ . The endpoints of each muscle remain fixed, while their lengths and orientations change dynamically. Forces and moments at the elbow joint were calculated using a dynamic multibody system formulation, in which the upper arm (humerus) remains stationary in a vertical configuration and only the forearm moves.

2.2. Calculation of Forces and Torque at the Elbow Joint

We perform an inverse dynamic analysis [27] to compute the unique joint forces and moments at the elbow for the desired curl, which will then be used to estimate the muscle force distribution. Although the joint moment has contributions solely from the muscles, the joint forces include contributions both from the muscles and from interactions between the bones.

The Newton-Euler equations governing the motion of the multi-body system are generally encapsulated by

$$M\ddot{q}=G(q,\dot{q}),$$

(1)

where, in our case,

$$\left[M\right]=\left[\begin{array}{ccc}m& 0& 0\\ 0& m& 0\\ 0& 0& \mu \end{array}\right]$$

(2)

is the inertia matrix (m is the mass of the forearm and $\mu$ is the moment of inertia specified in Table 1), and

$$\left\{q\right\}=\left(\begin{array}{c}x\\ y\\ \theta \end{array}\right)$$

(3)

is the vector of generalized coordinates. Here x and y are the coordinates of the center of mass (CM) of the forearm and $\theta$ is the angle that the body-fixed vector ${\widehat{a}}_1$ makes with the global vector ${\widehat{e}}_1$ as shown in Figure 1. Quantities $\dot{q}$ and $\ddot{q}$ are the generalized velocities and accelerations, respectively, and $G(q,\dot{q})$ represents the generalized forces that make up the equivalent force system at the center of mass, of the constraint (joint) forces and moment that we intend to compute, as well as the applied loads. In our case,

$$\left\{G\right\}=\left(\begin{array}{c}0\\ -mg-W\\ -W(L-L_{cg})sin\theta \end{array}\right)+\left(\begin{array}{c}{\lambda }_1\\ {\lambda }_2\\ T-{\lambda }_1{L}_{cg}cos\theta -{\lambda }_2{L}_{cg}sin\theta \end{array}\right),$$

(4)

where the first part captures the contribution of the applied loads (the load W carried by the arm and the weight $mg$ of the forearm itself) and in the second part, the Lagrange multipliers ${\lambda }_1$ and ${\lambda }_2$ are the forces at the elbow joint and T is the joint torque that we intend to compute. Here L is the length of the forearm and $L_{cg}$ is the distance of CM from the elbow specified in Table 1. The elbow, considered a revolute joint, imposes a set of constraints that are appended by a driving constraint of constant angular velocity $\omega =\dot{\theta }$ of the forearm. These constraints are described by the following set of expressions that are all set to zero.

$$\{\Phi \}=\left[\begin{array}{c}x-L_{cg}cos\theta \\ y+L_{cg}sin\theta \\ \theta -\omega t-{\theta }_{\circ }\end{array}\right],$$

(5)

where, ${\theta }_{\circ }$ is the starting elbow angle, which is ${\displaystyle \frac{\pi }{12}}$ and t is the time. Setting each of the expressions of $\{\Phi \}$ to zero, the generalized coordinates q that describe the arm configuration are solved at each time instant. Then, the generalized velocities $\dot{q}$ and generalized accelerations $\ddot{q}$ are solved at each time instant using

$$\dot{q}={\Phi }_q^{-1}{\Phi }_t$$

(6)

and

$$\ddot{q}={\Phi }_q^{-1}\gamma ,$$

(7)

where

$$\gamma ={\left({\Phi }_{q}q\right)}_q\dot{q}+2{\Phi }_{qt}\dot{q}+{\Phi }_{tt},$$

(8)

${\Phi }_q=Jacobian\left(\Phi ,q\right)$ is the constraint Jacobian and likewise ${\Phi }_t=Jacobian\left(\Phi ,t\right)$, ${\Phi }_{qt}=Jacobian\left({\Phi }_q,t\right)$, and so on depending on the subscripts in the notation.

The time series of q, $\dot{q}$ and $\ddot{q}$ solved using the above-mentioned steps are then substituted in Eq. (1) to get time series of the joint forces ${\lambda }_1$ and ${\lambda }_2$ and the torque T. To estimate the muscle forces that contribute to these joint forces and moment, the muscle model needs to be considered, which are presented in the next section.

2.3. Muscle-Force Models

In this study, we evaluated the performance of three skeletal muscle force generation models, listed below in order of increasing complexity:

• Simple Force Model [28]:
  This model considers muscle force as a vector directed from its insertion to its origin, without considering the effects of the force-length relationship of the skeletal muscle or the speed of muscle contraction.
• Hill-Type Active Model[20]:
  This model incorporates active muscle behavior influenced by muscle fiber length and contraction speed. [20]. The muscle force equation is:

  $$F_m\left(t\right)=F_{max}^{iso}·a_m\left(t\right)·f_a\left(l_m\right)·f_v\left(l_v\right),$$

  (9)

  where $F_{max}^{iso}$ represents maximum isometric force, $a_m\left(t\right)$ is muscle activation, $f_a\left(l_m\right)$ captures the active force-length relationship and $f_v\left(l_v\right)$ captures the force-velocity relationship. This model provides a dynamic view of muscle performance under varying physiological conditions.
• Hill-Type Active & Passive [29]:
  This comprehensive model includes both active and passive muscle behaviors accounting for muscle fiber length and muscle contraction speed [29]. The muscle force is given by:

  $$F_m\left(t\right)=F_{max}^{iso}.[a_m\left(t\right).f_a\left(l_m\right).f_v\left(l_v\right)+f_p\left(l_m\right)].$$

  (10)
  Here, the additional term $f_p\left(l_m\right)$ represents the passive force-length property, adding another layer of complexity by considering the intrinsic properties of muscle fibers.

The active and passive force-length properties are expressed in [20,30] as

$$f_a\left(l_m\right)=e^{\left(c{\left|{\displaystyle \frac{l_m-l_{opt,m}}{l_{opt,m}w}}\right|}^3\right)}$$

(11)

and

$$f_p\left(l_m\right)=e^{\left(10{\displaystyle \frac{l_m}{l_{opt,m}}}-15\right)}$$

(12)

respectively, where $l_m$ is the muscle fiber length, $l_{opt,m}$ is the optimum muscle fiber length to produce maximum force, w determines the width and is set as 0.4 and c is a coefficient set as $ln\left(0.05\right)$.

The force-velocity property for muscle shortening is given as [30]:

$$f_v\left(l_v\right)={\displaystyle \frac{-12l_{opt,m}-\dot{l_m}}{-12l_{opt,m}+K\dot{l_m}}},\dot{l_m}<0$$

(13)

where ${\dot{l}}_m$ is the muscle lengthening velocity and K is the curvature constant and is set to 5. Both Active only and Active & Passive Hill-Type muscle models rely on the force-length relationships of both muscle fibers as well as tendons. These relationships are challenging to model due to their nonlinear nature and dynamic changes during contraction. Accurately determining muscle fiber length at each time step is particularly challenging, as the total muscle length includes both muscle fibers and tendons, which change non-linearly. To address these issues, our approach uses known optimal muscle lengths and angle ${\alpha }_i$ as starting points to model muscle and tendon behavior.

Additionally, simplifications and assumptions regarding elasticity and stiffness are applied to represent muscle and tendon behaviors. These modifications simplify the model while retaining essential characteristics of the muscle-tendon system, thus allowing for the development of models that can reasonably estimate muscle forces and their interactions within the system, despite the complexities of nonlinear force-length relationships and variable lengths during contraction. The following adaptations are often made to address specific modeling challenges.

• Constant tendon length - In this approximation, muscle lengths exceed realistic bounds, and force-length properties are inaccurately represented at the beginning and end of the motion.
• Linear muscle contraction - In this adaptation, the muscle length is assumed to change at a constant rate. As a result, however, the force-velocity value is constant and is not correctly represented.
• Linear tendon length change - In this adaptation, while the rate of tendon length change remains constant, muscle length changes non-linearly, providing more realistic force-length and force-velocity values.
• Exponential tendon length change - In this one, the rate of change of tendon length varies exponentially, offering the most accurate representation of the muscle model.

A comparative analysis between the Linear and Exponential Tendon Length approximations reveals that the force-length and force-velocity values closely match the actual data observed during the bicep curl motion from ${15}^{\circ }$ to ${120}^{\circ }$ [31]. Consequently, we chose the Linear Tendon Length approximation for our simulation case study. This choice simplifies the modeling process while still producing results that closely align with the observed force-length and force-velocity characteristics during the curling motion, and thus effectively balances accuracy with computational efficiency, making it suitable for our intended purpose.

2.4. Muscle-Force Constraints and Cost Functions

The single DoF of elbow flexion represented by angle $\theta$ is controlled by three muscles (bic, bra, brd) which independently exert forces to produce torque around the elbow joint. This setup allows infinite combinations of muscle force distributions to achieve the same movement. This over-actuation presents a challenge in uniquely identifying the muscle forces that must undergo some constrained optimization. Although the net joint torque T required to perform the motion can be uniquely derived from the way $\theta$ evolves using the inverse dynamic approach, identifying the specific contributions to the torque of the forces exerted by each muscle group on the forearm (${\overrightarrow{F}}_{bic},{\overrightarrow{F}}_{bra},{\overrightarrow{F}}_{brd}$) is challenging due to the lack of sufficient equations. These forces are related to the torque T as follows:

$${\overrightarrow{r}}_{bic}\times {\overrightarrow{F}}_{bic}+{\overrightarrow{r}}_{bra}\times {\overrightarrow{F}}_{bra}+{\overrightarrow{r}}_{brd}\times {\overrightarrow{F}}_{brd}=T{\widehat{e}}_3,$$

(14)

where each term is a cross-product of position vector ${\overrightarrow{r}}_i$ of $i^{th}$ muscle’s insertion point on the forearm relative to the elbow joint with the corresponding muscle force. To solve the inverse problem of finding muscle forces from the above equation, optimization approaches and constraints that appropriately narrow down the range of possible force combinations are needed.

One constraint applied in this study involves the maximum isometric force that each muscle can produce, which sets upper limits for muscle forces. This constraint is expressed through the inequality:

$$0\le |{\overrightarrow{F}}_i|\le F_{max,i};i\in [bic,bra,brd],$$

(15)

where $F_{max,i}$ represents the maximum force that the $i^{th}$ muscle can exert. The values for the maximum isometric forces for the bic, bra, and brd muscles were taken from existing literature [29]. Depending on the muscle model used, additional constraints on muscle force combinations vary: (a) for the Simple Force Model, no further constraints are applied; (b) for the Active Hill-Type Model, muscle forces are constrained by active force-length and force-velocity relationships, as specified in Eq. (11,12); (c) for the Active & Passive Hill-Type Model, the constraints include active force-length, passive force-length, and force-velocity properties, detailed in Eq. (11,12,13).

These constraints restrict the range of possible solutions, and so a cost function must be applied to find a unique solution. This study evaluates the three different cost functions listed in Table 2: Cost function $J_1$ minimizes the total muscle force exerted [14,16,19], $J_2$ minimizes the total work done by the muscles [17], and $J_3$ minimizes the total stress developed in the muscles, the stress being the force per physiological cross-sectional area (PCSA) of each muscle[9,14,18,19]. The values for $PCS{A}_{m,i}$ were obtained from references [32,33], ensuring that the calculations are based on established biomechanical data.

• F_m,i: magnitude of force exerted by the ith muscle;
• Δl_m,i: change in length of the ith muscle;
• PCSA_m,i: physiological cross-sectional area of the ith muscle.

2.5. Model Validation

For this study, we conducted validation by comparing the calculated moment arms and muscle moments for each muscle with data from previously published studies. In particular, the comparison utilized moment arm data from a diverse subject pool (both male and female) captured during elbow flexion movements, starting from ${25}^{\circ }$ and extending to ${110}^{\circ }$ as documented in reference [34]. Additionally, muscle moment data spanning an extended range of elbow flexion from $0^{\circ }$ to ${120}^{\circ }$ were used for further validation [35]. Our model specifically simulates elbow flexion from ${15}^{\circ }$ to ${120}^{\circ }$, making it compatible with the ranges used in these studies for a comprehensive comparative analysis. The angular velocity of the elbow joint within our model simulation was held constant at $\omega =1.5\mathrm{rad}/\mathrm{s}$, which is slow enough to reduce dynamic effects and expect the results to be in the ballpark with the benchmarks cited above.

To get an effective dimensionless comparisons between the results, the moment arm and muscle moments were normalized.

3. Results and Discussion

This section details the results of our simulations, analyzing the generated joint torques and muscle forces during elbow flexion. We begin by exploring the torque dynamics required to rotate the elbow from ${15}^{\circ }\mathrm{to}{120}^{\circ }$ in Section 3.1, highlighting changes throughout the motion. Subsequently, in Section 3.2 we evaluate muscle forces under three modeling scenarios: the Simple Force Model, which serves as a baseline; the Hill-Type Active Model, which incorporates active force-length relationships for a more realistic muscle behavior; and the Hill-Type Active and Passive Model, which further integrates passive force properties to enhance model accuracy. In Section 3.3, we validate our model by comparing our predictions with published data, focusing on muscle moments and their moment arm. This comprehensive analysis not only confirms the efficacy of the employed models and cost functions but also highlights their practical implications in clinical and sports applications as discussed in Section 3.4. The findings underscore the trade-offs between model complexity and computational efficiency, emphasizing how even simpler models can achieve a high degree of accuracy when appropriately configured.

3.1. Joint Torque

Figure 2 illustrates the torque required to rotate the elbow with a constant angular velocity from ${15}^{\circ }$ of to ${120}^{\circ }$ of flexion over $1.22\mathrm{s}$. The torque increases with increasing elbow flexion until the arm reaches ${90}^{\circ }$, where the external load W becomes orthogonal to the forearm. The torque then decreases as the joint rotates past ${90}^{\circ }$ to its final position as expected.

3.2. Muscle Forces

This subsection evaluates muscle forces during elbow flexion in three muscle models, highlighting the impact of small changes in tendon length. We start with the Simple Force Model, assessing basic force generation, then advance to the Hill-Type Active Model which integrates active force-length relationships, highlighting how slight tendon variations affect muscle dynamics. Finally, the Hill-Type Active and Passive Model incorporates passive properties, demonstrating detailed effects of tendon adjustments on overall muscle force accuracy.

3.2.1. Simple Force Model

By disregarding the constraints on muscle force generation imposed by the physiological force-length relationship, the Simple Force Model assumes that each muscle is capable of producing as much force as is necessary to drive the joint throughout the full range of elbow flexion. Therefore, the Simple Force Model allows us to assess how the choice of cost function exclusively influences the muscle force distribution. In particular it reveals that the choice of cost function significantly influences the distribution of muscle forces in biomechanical simulations as described below.

The results of the Simple Force Model, as shown in Figure 3a, indicate that minimizing the Sum of Forces (cost function $J_1$), the brachioradialis (brd) produces the largest forces, followed by the biceps brachii (bic), and then the brachialis (bra). Figure 4 illustrates that these force magnitudes are directly related to their respective moment arms. This outcome is predictable because muscles with longer moment arms have a mechanical advantage in generating torque. Consequently, these muscles tend to bear a larger share of the load, as increasing force in muscles with shorter moment arms contributes less effectively to the overall torque, making them less efficient in terms of torque production.

The Sum of Work cost function ($J_2$) aims to minimize the total mechanical work done by the muscles. Given that work is the product of force and the muscle’s change in length, muscles that change length less will need to exert more force. This function shows how biomechanical efficiency can be modeled and optimized by balancing force output against changes in muscle length. Among the three muscle groups, bra is the shortest, and so is the change in its length during flexion, followed by bic and brd. Hence, as seen in Figure 3a, bra exerts the greatest force, with bic next, and brd exerting the least.

The Sum of Stress cost function ($J_3$) aims to minimize total stress, where the stress in each muscle is defined as the force exerted by it per unit of physiological cross-sectional area (PCSA) of the muscle. Muscles with larger PCSAs can therefore sustain higher forces. This cost function is particularly useful for understanding how structural muscle properties like PCSA influence the biomechanical performance and force capabilities of different muscles. Among the three muscle groups, bra has the largest PCSA ($5.6{\mathrm{cm}}^2$), closely followed by bic ($5.1{\mathrm{cm}}^2$), and brd ($1.2{\mathrm{cm}}^2$) [32,33]. Since the PCSA for bra and bic are comparable, the muscle with the longest moment arm, bic, is leveraged to exert the highest force, followed by bra, and brd, as shown in Figure 3a.

Overall, using the Simple Force Model to evaluate these cost functions illustrates how each function impacts muscle force predictions differently based on mechanical and physiological principles.

3.2.2. Hill-Type Active

The addition of the active force-length relationship of skeletal muscle to the model limits the maximum force that each muscle can exert during contraction. As seen in Figure 6, when accounting for the force-length relationship for each muscle, the bic, bra, and brd can produce maximum active force at approximately ${95}^{\circ }$, ${90}^{\circ }$, and ${50}^{\circ }$ of elbow flexion, respectively. As muscle length shortens (higher flexion angles) or lengthens (lower flexion angles) relative to this optimum length, the maximum possible force that each muscle can produce decreases accordingly (Figure 5). For small flexion angles, the limits on maximum force production are much more pronounced for the bic and bra muscles than for the brd.

Indeed, when compared to the Simple Force model, the forces produced at the start of each Hill-type active simulation are significantly lower for the bic and bra muscles for all three cost functions, despite both muscles exhibiting relatively high levels of activation. These effects are particularly pronounced for the Sum of Work and Sum of Stress cost functions, where the bic and bra muscles are at maximum activation as shown in Figure 3b. On the contrary, the brd forces are higher than predicted by the cost functions alone at these angles even while experiencing less than $20\%$ activation.

The effects of force-length relationships on muscle recruitment continue throughout the curling motion. As the flexion angle approaches the ${90}^{\circ }$ where the bic and bra are nearer their optimal length, their force production increases even as their activation level decreases. At the same time, the brd shortens past its optimal length so that its force-producing capabilities drastically decrease. As such, the brd rapidly approaches maximum activation, its force production levels, and eventually drops. The sudden decrease in brd force production coincides with rapid increases in bic and bra force production as they approach their optimal length, and the applied moment about the elbow reaches a maximum.

3.2.3. Hill-Type Active and Passive

The passive force-length relationship further tunes the initial results by including additional force potential in the fully extended arm position, similar to an extended spring, as shown in Figure 5. The addition of the passive force-length relationship of each muscle plays an important role in muscle forces and activation, especially as the stretch ratios moves beyond $1.1$. When the arm is in fully extended position, all three muscles are in a stretched configuration that generates large passive forces that exponentially decrease to low values when the muscles reach their optimum length. Due to this additional force the muscles can exert, the forces trend for bic and bra change for the passive model and the forces value are higher for sum of work and the sum of stress cost functions, where these muscle forces dominate shown in Figure 3c. Similarly, the decrease in activation of the bic and bra muscles at small flexion angles is partially due to the large passive muscle forces at these angles. Once each muscle has shortened past its optimum length, the force and muscle activation trends become similar to the active-only model, since passive force contributions become negligible. The inclusion of both active and passive components allows the model to capture the complex, nonlinear behavior of muscles, which became the deliberate choice for further analysis and validation.

3.3. Model Validation

The comparisons between the normalized moment arm data Figure 7 reveal a similar trend of change in moment arm for each muscle with respect to the angle of elbow flexion. The normalized moment arm for the three elbow flexors gradually increases to its maximum, around ${90}^{\circ }$ where the moment is at its peak, then decreases. We suspect that the slight variation arose from the variation of muscle length, PCSA value, muscle placement, and relative orientation among the individual subjects taken into account from the literature. In Figure 8c, the energy cost function $\left(J_1\right)$, the work cost function $\left(J_2\right)$ and the stress cost function $\left(J_3\right)$ provide normalized moments that align well with the observed moments for the brd muscle. For the bra muscle (Figure 8b), the work and stress cost functions show fewer deviations from the results in [35] compared to the force criterion. Regarding the estimation of muscle moment, the stress cost function for the bic muscle (Figure 8a) appears to show the closest agreement with the findings in [35].

It is worth noting that the stress cost function captures the increasing trend of the muscle moment with greater accuracy for all the muscles. These observations are consistent with the statistical analysis (Table 3) done on the plots with reference to the muscle moments plots obtained by [35] using Adam’s software. The stress cost function ($J_3$) has the least mean squared error for the bicep and brachoradialis, while the work cost function ($J_2$) has the least error for the brachialis muscle. The mean squared error of $J_3$ is the least, second to $J_2$ for Brachialis.

This characterization showed that all three simple cost functions are reasonably accurate in estimating muscle forces for brachioradialis, while for brachialis, the work and stress cost functions gave a fair estimate and the stress cost function provided a fitting estimate for the bicep. Overall, the combination of stress cost function along with the Hill-type muscle model with active and passive force length relationship provides the most reliable estimation of muscle forces for the three muscle groups during an elbow flexion. This simple cost function sensibly captures the increase in muscle force to its peak. The deviations observed beyond the peaks can be attributed to variations in the task configurations and the estimated muscle parameters for different subjects. The trend of increasing muscle force is an essential muscle behavior, and its estimation can be beneficial in the aforementioned domains where accurate force measurements are not the primary focus, such as performance evaluation, injury prevention, muscle activation assessment, etc.

3.4. Model Applications and Future Work

Our approach focuses on achieving a delicate balance between accuracy and usability, acknowledging that the pursuit of the utmost precision is not always paramount. This equilibrium drives us to explore scenarios where a reasonable level of accuracy suffices without necessitating the use of highly complex models. This balance or trade-off is useful in situations where incorporating EMG data could introduce complexity or when these data are simply not accessible. In essence, the approach does not rely on intricate EMG data, ensuring that it remains accessible and valuable in various scenarios. This exploration holds significant value in practical applications where a detailed biomechanical analysis may not be the primary focus. For example, a) Physical Therapy and Rehabilitation monitoring, where the primary goal may be to monitor patient progress and track changes in muscle activation rather than achieving precise muscle force estimation. b) Exploratory studies of biomechanical research, specifically while developing new methodologies for muscle activation patterns and trends exploration. c) Athletic Training and Performance Assessment, where the focus may be on performance assessment, movement optimization, and injury prevention, and d) Basic Physiology and Exercise Science). In future work, this study can be extended to formulate models for other body segments or joints and for different motions, such as knee extensions, abduction, adduction, rotation, etc.

4. Conclusions

This study demonstrates that the use of simple cost functions does not significantly compromise the accuracy of muscle force estimates, provided that they are combined with a sufficiently representative muscle model. The choice of the physiological muscle model is more important than the cost function as a more accurate muscle model captures the nonlinear behavior of the muscle, and the cost function provides an optimal strategy. Relatively simple mathematical cost functions can drive muscle force optimization and capture the essential aspects of muscle behavior, excluding the requirement for EMG data. This approach is highly beneficial in cases where the respective EMG data is simply not accessible. It can be effective and favorable for estimating the behavior of deep muscles, as EMG data for these muscles are difficult to obtain and prone to measurement errors. Furthermore, considering separately the cost function optimization and the muscle model allows for evaluation of their impact independently on the forces estimation as well as their combined effects.