Biceps brachial neural control simulations using OpenSim

The aim of the study was to evaluate the potential of OpenSim to test the neural control of the biceps brachii in flexion and extension movements. We used the OpenSim API and developed a code in Matlab to simulate the flexion/extension in the sagittal plane of an elbow with one kinematic degree of freedom. The neural command was then modified and executed with different parameters, over 3 simulations. To simulate the elbow flexion/extension, two bodies were created, to represent the humerus and radius bones, and an actuator, to represent the biceps muscle. The parameters changed at each simulation were: the times of the minimum/maximum excitation of the biceps brachii muscle and the minimum/maximum values of excitation. The values chosen for simulation 1 of maximum and minimum excitation time were 0.5s and 3s, while the excitation values were 0.3 and 1, whereas in simulation 2 the chosen time values were 3s and 7s, with excitation values equal to first simulation. Finally, simulation 3 had as time parameters 1s and 10s and the maximum and minimum excitation values being 0.5 and 1. In simulation 1, a smaller extension angle value was observed when compared to simulation 2, but with similar force values in both simulations. In simulation 1, it is possible to observe that in a shorter excitation time, a muscle contraction response occurs, interpreted by the resulting force. In simulation 3, the same strength level was observed when compared to simulation 2, however, for a longer period of time. We conclude that with the 3 simulations performed, it was possible to analyze the behaviors of biceps force and elbow joint angle according to the neural stimulus simulated by the muscle excitation values. In addition, simulation 1 can be used in future studies that seek to represent a neurological patient with


I. INTRODUCTION
Computational models have been used as tools in studies that address the neuromotor system, in part, because they allow the measurement of values that are difficult to measure in vivo, such as: tendon stretch tension, individual muscle activation, joint moments, among others [1,2]. In addition, computational models allow for a reduction in the excessive repetition of data collection with volunteers undergoing the study. Thus, favoring a more accurate data collection and reducing the risk of injuries to patients.
OpenSim is an open source software that makes it possible to model the biomechanics of movement in a realistic way [3] and, therefore, has several applications in the clinical or sports area. One such application is the investigation of neural control mechanisms. This has already been done by other authors using a musculoskeletal model containing muscle spindle, implemented in OpenSim [4]. The muscle spindle is the sensor responsible for providing the spinal cord with information on muscle stretch length and speed to produce stretch reflexes that also depend on commands or modulations sent by the brain. It is noteworthy that these reflexes cause changes in muscle activation [5].
In healthy condition, proprioceptive information is processed by the nervous system to correct errors under unpredictable environmental conditions, mediated primarily by feedback control [for a review, see 6]. Feedback control for instantaneous corrections was proposed to be implemented through changes in stretch reflex thresholds, regulated by the corticospinal system (CS) [7,8]. This guarantees positional and velocity stabilization despite of unexpected external perturbations. In addition, under consistent and predictable environmental conditions, proprioceptive information is processed to ensure trajectory accuracy and movement planning during motor tasks [6]. According to the studies mentioned above, it would be expected that a brain injury, such as a stroke, would produce sensorimotor deficits in both upper limbs during motor tasks.
In the presence of post-stroke people, the CS's ability to regulate muscle activation through changes in stretch reflex thresholds is impaired and contributes to sensorimotor impairments such as spasticity [7].
In the literature, a model has been proposed to represent the reflex pathways of the elbow. For example, Hidler and Schmit [9] aimed to understand the relationship between the speed of elbow flexion with the magnitude of the active reflex force. The authors hypothesized that the observed force plateau could be explained by an inhibitory force feedback pathway. The study compared simulations of elbow reflex pathway models with data collected from patients with some level of post-stroke spasticity. The authors used a model that contained two separate feedback pathways, one representing the monosynaptic stretching reflex originating from muscle spindle excitation and the other representing the inhibition of force feedback originating from force-sensitive receptors. It was found that the inhibition of force feedback altered the response of the stretch reflex, resulting in a force response that followed a sigmoidal shape, occurring in greater amplitudes of the spastic elbow [9].
The present study aimed at evaluating the viability of using Matlab-OpenSim interface to generate muscle forces results comparable to physiological behavior, by implementing a possible mechanism of neural control of the brachii biceps in flexion movements, from commands generated based on the muscle model proposed by Millard, et al. [10]. The model proposed for our investigation uses a parameter of muscle response, providing more reliable values of strength in a real situation. In this sense, force feedback relationships are extremely important to better understand the pathological mechanisms of spasticity [9].

II. METHOD
All programming was performed in MATLAB®, 2019b, (MathWorks Inc, Natick, Massachusetts) using the Matlab-OpenSim interface tools and the OpenSim Moco toolkit available in OpenSim 4.0. The OpenSim Moco, is a software toolkit that allows one to solve optimal control problems using the musculoskeletal models defined with OpenSim. The Opensim and OpenSim Moco functionalities are accessed through Matlab as a Java Class Library.
We, thus, created a simplified model of the upper limb to simulate the elbow flexion-extension, starting from a OpenSim codding tutorial, available at opensim-core git hub (https://github.com/opensim-org/opensim-core), and modifying it to represent the human motion more realistically. The model available in the code included two rigid segments, representing the upper-arm and the forearm, which were connected by a pin-joint, allowing for one rotational degree of freedom (DOF), that could account for the elbow flexion-extension in the sagital plane. The model was actuated by one hill-type muscle, that could be used to represent the biceps Brachii, an elbow flexor. Muscle force generation and activation dynamics were described according to Millard, et al. [10].
A detailed description about the model modification as well as the simulations performed are described in the following sections.

A. Biomechanical Architecture
Within the example-code the segments were already created with the class "Body" and the connection between them was made using the class "PinJoint". The segment dimensions and inertia properties, however, were not properly adjusted to represent an average adult. Thus we properly scaled the inertial properties and segment dimensions of both segments to represent a 60kg weight and 1.70m height person. The inertial characteristic, namely, upper and forearm mass an moment of inertia, were estimated from Winter [11] model and the respective segment longitudinal length were estimated according to Winter [11] model (see Table 1).

A. Muscle Model
Muscle Model was set using the class "Millard2012EquilibriumMuscle", and have the main equation (1) of muscle dynamics: Where 0 , is the maximal isometrical force, is muscle activation which maximum possible value is 1, (̃) is the force-length relationship, (̃) is the force-velocity relationship, (̃) is the force-length relationship of passive structures, is the pennation angles of muscle fibers, (̃) is the tendon force-length curve and is the coefficient of dumping fibers [10]. Maximal Isometric Force was set as 109 N, optimal fiber length 0.34 m, Tendon slack length 0.125 m and pennation angle 0 rad. All parameter values were chosen or calculated using the equations and tables of Giat, et al. [12] for the biceps brachii.

B. Controller
Activation (a) is described as a percentage of muscle fiber recruited when a neurological input (u) comes from the upper neural centers. Equations that represent this relation are described below: Controller settings come from the class "PrescribedController" which allows applying three configurations: tendon properties, type of signal, and dumping fiber term configurations. Thus, we assume respectively an elastic tendon, a system of dynamic activation, and the default damping influence at the muscle fibers, which is 0.1. Activation and deactivation time constants are 0.01 and 0.04, respectively. The step function has been chosen to describe the neural input, which is a parameter required by the "prescribeControlForActuator" class.

C. Simulations
Three simulations were performed, with changes in the initial and final excitation times, and , respectively; excitations maximum and minimum values   Table 1. III. RESULTS Fig. 2 shows the behavior of excitation, elbow joint angle and force produced from the three simulations simultaneously for better comparison. In the 0.5s is registered the maximum angle during simulation (74º), to then stop in the position of as 66º in the end of the simulation. Regarding these changes, the outputs describe the muscle behavior and activation dynamics.     In the first second of the simulation number 1 (Fig 4), we can correlate the minimum excitation period with the huge amplitudes of angle variation and force produced. Assuming that correlation, the curves of force production and angle variation shows a physiological behavior when force increase, and the angle assumes minimum valuers of flexion. Fig 3-B, as in the simulation one shows the history of movement that the model does in 10s. However, three frames were simulating, with the same initial condition (90º) the movements start until achieve a lower angle (67º) in the 0.3s and then increases to 74° (0.5s) to finally, like in simulation 1, stop in the position of 66º.

B. Simulations number two
In

C. Simulations number three
The trajectory of the movement in simulation number 3 is shown in Fig 3-C. With three frames, the initial condition was set (90º) and close to the 0.2 s of simulation, the minimum valuer of angle is achieved (59º). Then the angle increases to 73º to decrease again, and the arm achieve the final position (66°).
For the last simulation a higher time of rise was set ( Table  2). Minimal values of angle were registered for higher excitation compared to the other simulations, similarly to simulation 2 correlating the force and short angles. The third simulation achieved the highest and lowest force values of the three simulations (16,7 N and -1.041x10^-16 N).

IV. DISCUSSION
Hill models are often used in simulations. Even though they represent a simplification of the physiological properties of a muscle contraction, these models effectively meet the demands of the simulation [2]. These models are essentially used when they involve many muscles, which then have their calculations simplified, when compared to models that try to contemplate the entire physiology of muscle contraction [10].
According to our results, at the end of the 3 simulations (Fig. 4, 5, 6) it was possible to observe that changing the biceps muscle excitation time values, as well as its minimum and maximum value through a step function, the Muscle responds to this stimulus by changing the angle of the elbow joint and muscle fiber force. Thus, increasing the time for the minimum excitation value and the time for the maximum excitation value from simulation 1 to 2 ( Fig.  4 and 5), with the same maximum and minimum excitation values, it was noticed that simulation 1 resulted in a lower extension angle value when compared to simulation 2, but with similar muscle fiber strength values in both simulations. Furthermore, in simulation 1 and 2 ( Fig. 4 and 5) they are similar to represent the force generated by the biceps, showing that they reach a maximum force value of 9.4 N in 0.1 s, whereas in simulation 3 it has a maximum value of 16.7 N in 0.4 s. These results suggest that using the parameters of simulation 3 it is possible to generate the same level of muscle strength when compared to simulation 2, for a longer period of time.
The results generated by the three simulations in this study suggest that the excitation parameters used in the three simulations were considered satisfactory to physiologically represent the muscle behavior showing the relationship strength x length [13]. The strength x length relationship proposes that strength is dependent on the length of the muscle, where studies show that this behavior is related to the two main proteins of the contractile muscle element: actin and myosin. These proteins are present in the sarcomere, and overlap and slide with each other, based on the theory of cross-bridge contraction [14]. Simulations with changes in time parameters and joint angles are useful to understand the neural mechanism of excitation and contraction in health and clinical populations.
The dynamics of muscle tone control originates from afferent pathways (proprioceptive feedbacks coming from the neuromuscular spindle) and efferent pathways (higher command centers via tonic stretch reflex). However, in individuals affected by stroke, this dynamic becomes altered as the higher centers that modulate the tonic stretch reflex (TSRT) are affected [15]. Such changes in TSRT modulation are directly related to the extension and affected areas that result in degrees of severity. In the case of post-stroke, this exaggerated response is closely related to the speed of force generation, resulting in spasticity [for a review, see 16]. In this sense, simulation 1 could be used in future studies that seek to represent a neurological patient with spasticity. Spasticity is characterized by muscle hyperactivity due to a velocity-dependent increase in tonic stretch reflexes that results from abnormal spinal processing of proprioceptive stimuli [17]. The spastic muscle responds to the stimulus much faster than a muscle without spasticity. In simulation 1, it is possible to observe that in a shorter excitation time there is a muscular contraction response, interpreted by the resulting force. Future studies are needed to further clarify the best parameters to represent a spastic muscle.
We conclude that with the 3 simulations performed, it was possible to analyze the behaviors of biceps force and elbow joint angle according to the neural stimulus simulated by the muscle excitation values.

Final considerations
In this research, a step function was used to describe the actuator, defining the way biceps excitation should behave over time, but it is possible to use other functions to define this stimulus and that are available in the OpenSim API documentation, such as the functions "PiecewiseConstantFunction", "GCVSpline" and "SimmSpline", which can be used in future research, analyzing the parameters needed to use them, how their curves behave, to later verify what this entails in the simulation behavior.