<figure><img src="image_1.png" /></figure>

HAL
open science

Modeling of Isometric Muscle Properties via
Controllable Nonlinear Spring, and Hybrid Model of
Proprioceptive Receptors

Mario Spirito

► To cite this version:

Mario Spirito. Modeling of Isometric Muscle Properties via Controllable Nonlinear Spring, and Hybrid Model of Proprioceptive Receptors. 2023. hal-04088329v2

HAL Id: hal-04088329

https://hal.science/hal-04088329v2

Preprint submitted on 31 May 2023

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Modeling of Isometric Muscle Properties via Controllable Nonlinear Spring, and Hybrid Model of Proprioceptive Receptors

Mario Spiritoª,

ª46 Boulevard Niels Bohr (LAGEPP - Bâtiment C.P.E.) - Université Claude Bernard Lyon-1, Villeurbanne, France

## Abstract

This work presents the macroscopic behavior of skeletal muscles seen from a system theoretic standpoint. Then, according to the data available in the literature, we propose a first Evaluation model for the muscle isometric (constant muscle length) force generation depending on the muscle length and neuronal excitation frequency. Such a model is instrumental in providing a more physics-inspired model of such an isometric force by exploiting a nonlinear spring description with controllable characteristics, such as stiffness and rest length. To conclude, we propose a hybrid dynamical filter model to describe the components of the sensory system. This is in charge of feeding back to the Central Nervous System the measurements of muscle length and its rate of change.

**Keywords:** Skeletal muscles isometric force modeling, Nonlinear variable spring, Hybrid dynamical modeling, Spindle-Ia & -II, Oculomotor system.

## 1. Introduction

Skeletal muscles, or simply muscles, are organs of the muscular system whose cells have the characteristic to produce force and movements when electrically stimulated. They are mostly attached to the bones of the skeleton by tendons, and thanks to their contractile and elastic properties they allow the relative motion of the bones of the involved articulations.

Each skeletal muscle consists of thousands of muscle fibers wrapped together by connective tissue sheaths, see any human anatomy book, e.g. [1][Ch.9]. The individual bundles of muscle fibers in a skeletal muscle are known as fasciculi. Each muscle fiber is comprised of several myofibrils containing many myofilaments. When bundled together, all the myofibrils get arranged in a unique striated pattern forming sarcomeres which are the fundamental contractile unit of a skeletal muscle. The two most significant myofilaments are actin and myosin filaments which are arranged distinctively to form various bands on the skeletal muscle [2][Ch.1], [3] and [4].

Skeletal muscles show up with neuronal innervations from the sensory-motor system. A set of neurons, the α-motoneurons, constitutes the neuromuscular junction. They are responsible, via their excitation frequency here labeled as α (measured in pulses per second), for the sarcomeres contractions, and hence they allow the generation of the muscle active force. On the other hand, another set of neurons, constituting the proprioceptive receptors, carries sensorial information from the muscle to the Central Nervous System. In particular, the transferred signals are the values of muscle length and rate of change (Spindles -Ia and -II), and the measurement of the force value at the

tendon level (Ib-Spindle or Golgi Tendon Organs (GTO)). We can see the muscle structure, from a system theoretic standpoint, as a multi-input multi-output system with the actuator (α-motoneurons), the plant (the muscles cells, whose contraction produces the exerted force on the joint), and the sensors (Spindles -Ia and -II, and GTO), and the CNS as a feedback controller. A schematic of the involved neurological connections is available in [5].

<figure><img src="image_3.png" /><figcaption>Figure 1: Muscle and Central Nervous System schematics of the connections. Picture taken from [5].</figcaption></figure>

Email address: spirito.mario@gmail.com (Mario Spirito)

Preprint submitted to International Journal of Robust and Nonlinear Control
The modeling of how muscles exert force started with the work of [6] in which the muscle description via a spring-like behavior is presented. In particular, the relationship takes into account a contractile element that produces muscle force and it is put in series to a passive elastic element representing the tendon model between the muscle fibers and the connected bones or tissues. This model has then been exploited to generate more detailed models [7] and [8] in which the authors describe the microscopic behavior of contraction behavior generated by the effects of sarcomeres and cross-bridges. In [9], Hill's model has been enriched with a 'damping element'. In [10], [11], and [12], experimental data has been exploited to construct an analytic expression of the total force generated by the muscles. In these works, the total force has been split into a length-dependent force $F_L$ and a velocity-dependent force $F_V$, both depending on the $\alpha$-motoneuron excitation frequency. The $F_L$ term is then split, similarly to Hill's model, into an active plus passive term. The active force describes the part of fibers contraction when excited and produces movements, while the passive force only refers to the elastic mechanical behavior of the muscle fibers.

Experimental data and models of isometric (fixed muscle length) and isokinetic (fixed muscle velocity) forces are available and presented in [10], [11], [12], [13], and [14]. A system theoretic perspective of the neuromuscular system has been presented in [9, 15, 16] in which they give a control-oriented interpretation of the involved parts. In [17], instead, they consider the optimal co-activation problem while minimizing the mechanical impedance of an antagonist-actuated biological joint (e.g. an elbow).

In recent years, with the development of the contraction theory in the control field, the contractile properties of the neuromuscular system (involving Hill's model equations) have been studied in [18]. In [19] an analytic model of complete and incomplete tetanus contraction has been proposed, and, according to the authors' knowledge, it gives the closest to-reality dynamical model of the muscle behavior.

The sensory part involved in the muscle system is related to the measurement of the plant variables, such as force at the tendon level (GTO), and muscle/fibers length and velocity (Spindles - Ia and -II). They have been studied in depth [5] and [20] and the many works have been collected in [21], see [22] for a more recent review on the topic. From a system engineering perspective, the GTOs dynamical model has been proposed in [23] by exploiting a proper transfer function between the force exerted on the tendons and the GTO neurons' firing frequency. For Spindles -Ia and -II, in [24], a static model has been proposed to describe the relationship among the muscle/fibers length, velocity, and the $\gamma$-motoneurons excitation to the spindles firing rate. Such $\gamma$-motoneurons are in charge of modulating the firing responses of the connected Spindle system.

To the best of the authors' knowledge, there is no dynamic model of the Spindles organs available in the literature.

The sensory part plays an important role in the movement control architecture. In particular, the presence of internal loops between the muscles and the spinal cord allows the possibility to actuate the reflexes [16], [25], and [26], such as stretch re-

flexes, reciprocal inhibition, and others. Finally, another relative element is the spinal cord. It connects the CNS and $\alpha$-motoneurons via the interneurons and it is unavoidably involved in movement control [27, 28, 29].

The paper is organized as follows. We conclude the introduction with some motivations for the work, then in Section 2 we give some authors' considerations about the available data, and in Section 3 we describe models available in the literature that has been obtained exploiting these data. While in Section 4 we describe the developed models for the isometric muscle force. In Section 5 we provide a hybrid modeling of the sensory system available in the muscles. We then provide an application of muscle modeling here proposed to the description of the oculomotor system, in Section 6 and describe the challenges related to it. Section 7 concludes the work by summarizing the results and describing possible future works and advances.

### 1.1. Motivations

This work aims to describe neuromuscular behavior in a control theory paradigm. In particular, we would like to better understand the purposes of the neuromuscular loops and how they are exploited in controlling limbs movements, see for example [30] for the application of a stabilizing PID action on a single-link biomechanical model or [31] to see how this available measurement feedback can be exploited in position and movement control. The reader can furthermore see [18] to have an insight into how waveform signals can be used to control delayed closed-loop systems, thus mimicking the spiking behavior of biological nervous systems.

We strongly believe that this modeling idea would be a good starting point to eventually deeply comprehend the purposes of the brain parts involved in learning and control of dexterous movements, see for example [32], [33], [34], [35], [36], and [37], and thus, e.g., describe and simulate pathologies effects affecting movements control.

The final objective will then be to have a complete picture and understanding of the structure used by the brain in controlling human dexterous movements and thus mimic such a strategy in the field of nonlinear control systems such as nonlinear output regulation as in [38]. The study interleaves with the concept of inverse and forward models introduced in neuroscience in the works [39], [40]. See also [41] for another application of these models. Thus an intermediate step of the study will be to better understand how these models are necessarily formed, trained, and updated in the brain and how are they involved in the control strategy.

## 2. Preliminaries on available experimental data in literature

In the muscle force, we can distinguish among three main components: a passive force $F_p$ that is always present due to the natural behavior of tissues and only depends on the muscle length, then we have the active static (isometric) force $F_L$ only depending on the muscle contraction (length $L$) and an active dynamic (eccentric and concentric) force $F_V$ also depending on
the rate of change of muscle contraction $\dot{L}$. The last two components are active terms and indeed they depend on the activation frequency $\alpha$ of the $\alpha$-motoneuron.

Data of isometric force available in the literature are mainly those in the figures of [11], [12], and [17], reported in this section for completeness. It is worth noticing that such data only

<figure><img src="image_4.png" /><figcaption>Figure 2: Normalized Isometric Force. Relationship between the static Force and the muscle length at different activation frequencies. Plot taken from [12].</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 3: Isometric Force. Relationship between the static Force and the muscle length at different activation frequencies. This plot is half useful due to the unknown length, but we can notice that there are different rest lengths (length at which the force is null) depending on the activation frequency. Plot taken from [17].</figcaption></figure>

span the activation frequency in the interval [0, 120] pps (pulse per second), because 120 pps is the maximum excitation frequency the $\alpha$-motoneuron can ‘apply’ on the muscle fiber.

At $\alpha = 120$pps the active force $F_L$ (i.e. the total force without passive $F_p$ and dynamics $F_V$ terms) has a maximum value indexed in literature by the symbol $F_o$. This maximum value $F_o$ is associated with a specific length value $L_o$ (which is not the maximum length value, but the length $L$ corresponding to the maximum force value $F_o$ at the maximum excitation frequency $\alpha = 120$pps). Almost all available data and proposed models in the literature are normalized with respect to these two values. Indeed we find force and length data of the form $F/F_o$ and $L/L_o$, respectively. Due to this fact, we also consider a normalized model for such normalized values.

A crucial point in our model development is the presence in the force characteristic of rest length values (i.e. values of $L$ for which the exerted force is null) that depends on excitation rates

<figure><img src="image_6.png" /><figcaption>Figure 4: Normalized Isometric Force model with correction term given by a negative passive force in [11]. Here the activation frequency is normalized with respect to its maximum value $\alpha = 120$pps</figcaption></figure>

$\alpha$, see Fig.3. And this fact leads us to consider a rest length function $L_o$ (which is not $L_o$ as above) depending only on the excitation frequency $\alpha$ and this might be used to describe the contracting behavior of the muscle. A major problem is that most of the data available in the literature lack these rest-length data which is fundamental for our model development. We thus exploit the same technique introduced in [11]. In practice, they first consider a model without the presence of rest length (based on Gaussian terms, as it will be clearer at the end of section 3), and then they propose to add a negative passive term to compensate for the lack of rest lengths values in the model and obtain the rest lengths points saturating the force at zero. The approach is very clear in Fig. 4.

In the next section, we provide an overview of the two models available in the literature.

### 3. Literature models comparison

In our opinion, the most relevant and complete models in the literature are those available in [12] and [11]. Here we give a short description of such models for the sake of completeness.

#### 3.1. Brown's et al model

In [12], as the Hill's model [6], the authors consider both the active isometric and the dynamic forces generators gathered inside a *contractile element* $F_{CE}$ along which passive term $F_{PE}$ describes the passive force (previously indicated as $F_p$, but we decided to stay faithful to the original nomenclature). The resulting total force $F_{TOT}$ presents as

$$F_{TOT} = F_{CE} + F_{PE} \quad (1)$$

and $F_{CE}$ is expressed as

$$F_{CE}(\alpha, L, \dot{L}) = R(\alpha)A(\alpha, L)F_L(L)F_V(L, \dot{L}) \quad (2)$$

where $R$ is a fiber recruitment percentage and describes the amount of fiber involved in the force generation depending on the activation frequency, $A$ is the term relating the muscle activation to the excitation frequency $\alpha$ and it would have values between 0 and 1 (the same values range has been considered
for R), $F_L$ and $F_V$ describe respectively the force-length relationship and the force-velocity, and their product provides the force expression of the muscle depending on the current length L, the velocity of contraction $\dot{L}$, and excitation frequency $\alpha$.

In their work, a simplified version of the force function is provided fixing $R(\alpha) = 1$ and making assumptions of quasi-statics (close to isometric) conditions during motion and of a starting length equal to the rest length, hence they move to a force expression of the type:

$$F_{CE}(\alpha, L, \dot{L}) = A_F(\alpha, L) (F_L(L)F_V(L, \dot{L}) + F_p(L)) \quad (3)$$

where

$$ 
\begin{aligned}
A_f(\alpha, L) &= 1 - \exp\left(-\left(\frac{\alpha}{0.56N_f(L)}\right)\right) \\
N_f(L) &= 2.11 + 4.16\left(\frac{1}{L} - 1\right) \\
F_L(L) &= \exp\left(-\left|\frac{L^{1.93} - 1}{1.03}\right|^{1.87}\right) \quad (4) \\
F_p(L) &= -0.02e^{13.8-18.7L} \\
F_V(L, \dot{L}) &= \begin{cases} \frac{-5.72 - \dot{L}}{-5.72 + (1.38 + 2.09L)\dot{L}}, & \dot{L} \le 0 \\ \frac{2.5 + 4.21L + 2.67L^2}{0.62 + L}, & \dot{L} > 0 \end{cases}
\end{aligned}
$$

In this model, one can notice that the length $F_L$ and velocity $F_V$ dependent force terms are multiplied one each other. This might lead to the interpretation that the velocity terms $F_V$ are only acting as scaling functions on the 'supporting' isometric force $F_L$. We rather prefer to split the isometric $F_L$ and dynamics $F_V$ terms to analyze their contribution individually in a spring-damper fashion. Moreover, notice that when $\dot{L} \to 0^+$ in the equation we have $F_V \neq 1$, thus the model loses the continuity with the static isometric properties of the muscle after a lengthening phase.

## 3.2. Winters' Model

In [11], Winters's work goes deep inside the basic muscle elements. Indeed in its work, we can find a microscopic description of the muscle components such as intrafusal and extrafusal elements with their relative damping and stiffness terms, and these are connected in parallel and series fashion. Of their work, we only focus on the isometric force that has been modeled as

$$F_{ISO}(\alpha, L_{CE}) = F_L(\alpha, L_{CE}) + F_p(L_{CE})$$

where, for the isometric force model, they exploit a Gaussian function with variable mean value

$$F_L(\alpha, L_{CE}) = \alpha e^{-\left(\frac{L_{CE} - 1.05 - (1-\alpha)0.2}{0.4}\right)^2} \quad (5)$$

where $L_{m0}$, in Winter's notation, is the optimal muscle length value $L_o$ (in this work notation), $L_{CE} = L_{mt} - L_{t0} - x_{SE}(L_{m0} + L_{t0})$ and $x_{SE}$ is the (dimensionless) extension of the extrafusal series of spring-element, $L_{mt}$ is the overall musculotendinous

length and $L_{t0}$ is the tendon rest length. The value of $\alpha$ is normalized with respect to its maximum value. The passive force is instead given by the expression

$$F_{passive} = \begin{cases} \frac{e^{\frac{3}{0.6}\left(\frac{L_{CE}}{L_{m0}} - 1\right)}}{e^{3-1} - 1}, & \frac{L_{CE}}{L_{m0}} \ge 1 \\ -\frac{e^{\frac{6}{0.7}\left(1 - \frac{L_{CE}}{L_{m0}}\right)}}{e^{6-1} - 1}, & \frac{L_{CE}}{L_{m0}} < 1 \end{cases} \quad (6)$$

Notice that, as touched above, the passive force has a negative term for $\frac{L_{CE}}{L_{m0}} < 1$.

We strongly believe that, for our study purposes, the functions exploited in the models presented in this section are very complicated to manipulate and use in control analysis. Moreover, the meaning of the parameters is very difficult to interpret from a physical point of view.

## 4. Proposed Models

Our main objective is to describe the isometric force ($F_L + F_p$) as a nonlinear controlled spring with variable stiffness $K$ and rest length $L_0$

$$F_{ISO}(\alpha, L) = K(\alpha, L)(L - L_0(\alpha)) \quad (7)$$

where $\alpha$ is the $\alpha$-motoneuron excitation frequency and $L$ is the muscle length. In order to obtain such a final model, we first designed an intermediate/evaluation model following Winters' considerations. Indeed, we exploited the negative passive term introduced in [11] to obtain the rest length values for the different excitation frequencies. This step allows us to then construct and fit the parameters of the final/control model (7).

### 4.1. Evaluation Model

In order to construct the intermediate/evaluation model we consider the isometric force as divided in active $F_L$ and passive $F_p$ terms, as exploited in [11]:

$$F_{ISO}(\alpha, L) = F_L(\alpha, L) + F_p(L).$$

Inspired, then, by the modeling approach in [12] we describe the active force curves as Gaussian functions with mean value $\mu$ and standard deviation $\sigma$ and amplitude $A$ as a function of the excitation frequency $\alpha$, as shown in (4.1).

$$ 
\begin{aligned}
F_{\text{active}} &= A e^{-\left(\frac{L-\mu}{\sigma}\right)^2} \\
A &= \frac{c_{11}\alpha^2 + c_{12}\alpha + c_{13}}{\alpha^2 + c_{14}\alpha + c_{15}} \\
\mu &= c_{21}e^{c_{22}\alpha} + c_{23}e^{c_{24}\alpha} \\
\sigma &= c_{31} + c_{32}\cos(c_{34}\alpha) + c_{33}\sin(c_{34}\alpha)
\end{aligned}
$$

By running a curve fitting algorithm we are able to fit the active isometric force with a very good order of approximation, i.e. the standard Euclidean norm of the residual vector is $|r|_2 = 0.051$, where $r$ is the stack of the differences between
<figure><img src="image_7.png" /><figcaption>Figure 5: Data comparison between the available data in [12] and the evaluation/intermediate model. r indicates the residuals vector.</figcaption></figure>

<table><caption>Table 1: Parameters of Evaluation Model active force.</caption><thead><tr><th><i>i</i></th><th><i>c</i><sub>i1</sub></th><th><i>c</i><sub>i2</sub></th><th><i>c</i><sub>i3</sub></th><th><i>c</i><sub>i4</sub></th><th><i>c</i><sub>i5</sub></th></tr></thead><tbody><tr><td>1</td><td>1.134</td><td>0.723</td><td>0.0504</td><td>11.3297</td><td>643.259</td></tr><tr><td>2</td><td>0.574</td><td>-0.0211</td><td>0.8477</td><td>0.0008</td><td>n.a.</td></tr><tr><td>3</td><td>0.4518</td><td>-0.01704</td><td>-0.1235</td><td>0.04898</td><td>n.a.</td></tr></tbody></table>

the available force data and the related value of the proposed model. In Fig.5 we plotted the active force comparison between the available data from [12] and the intermediate model estimation. The identified coefficients are reported in Tab.1.

Exploiting the data available in [11] and their negative passive terms considerations we propose a unified continuous model of the passive force term exploiting a sigmoidal and exponential function

$$F_p(L) = f_1 + \frac{f_2 - f_1}{1 + 10^{f_3(f_4 - L)}} + f_5 e^{f_6 L} \quad (8)$$

were the coefficient values $f_i$, for $i = 0, ..., 6$, can be found in Tab.2. The sigmoidal behavior in the passive term, see Fig.9,

<table><caption>Table 2: Parameters of Evaluation Model passive force.</caption><thead><tr><th><i>f</i><sub>1</sub></th><th><i>f</i><sub>2</sub></th><th><i>f</i><sub>3</sub></th><th><i>f</i><sub>4</sub></th><th><i>f</i><sub>5</sub></th><th><i>f</i><sub>6</sub></th></tr></thead><tbody><tr><td>0.00381</td><td>2.16582</td><td>8.16634</td><td>1.45388</td><td>-7.32</td><td>-7.4</td></tr></tbody></table>

does not have an experimental data correspondence but provides a plausible description of the physical limit of every elastic material, i.e. beyond the ‘boundary’ length value the solicited fibers start to break. Then the total isometric force will result from the saturation at 0 of all negative terms, as in

$$F_{ISO}(\alpha, L) = \begin{cases} F_{passive}(L) + F_L(\alpha, L) \\ 0, & F_{ISO} \le 0 \end{cases} \quad (9)$$

## 4.2. Control Model

Due to the high complexity of the intermediate/evaluation model (9), and in order to make the force relationship more in-

<figure><img src="image_8.png" /><figcaption>Figure 6: Passive force comparisons among the available data in green, the Evaluation model in blue, and the control model in magenta. The residual values for the Evaluation model has norm $|r_{EM}|_2 = 0.0077$ and for the Control model the norm is $|r_{CM}|_2 = 0.1706$</figcaption></figure>

tuitive from a physical standpoint, we exploited the intermediate model to find the rest length points and then the equivalent stiffness values. These data were fundamental to then fit the rest length $L_0(\alpha)$ and stiffness $K(\alpha, L)$ functions' parameters in terms of the $\alpha$ and $L$, to construct the final/control model. Thus, the isometric force function is given by (7), reported here for completeness,

$$F_{ISO}(\alpha, L) = K(\alpha, L)(L - L_0(\alpha)) \quad (10)$$

where both $L_0$ and $K$ are taken as polynomial functions in $\alpha$ and $L$

$$\begin{align} L_0(\alpha) &= \ell_1 \alpha^2 + \ell_2 \alpha + \ell_3 \tag{11} \\ K(\alpha, L) &= k_1(L) \alpha^3 + k_2(L) \alpha^2 + k_3(L) \alpha + k_4(L) \nonumber \end{align}$$

where the $k_i(\cdot)$, $i = 1...4$, a polynomial functions in $\alpha$

$$\begin{align} k_i(\alpha) = k_{i1}\alpha^6 &+ k_{i2}\alpha^5 + k_{i3}\alpha^4 + k_{i4}\alpha^3 + k_{i5}\alpha^2 \nonumber \\ &+ k_{i6}\alpha + k_{i7}. \tag{12} \end{align}$$

The relative coefficient values are reported in Tables 3 and 4, respectively.

<table><caption>Table 3: Parameters of Control Model rest length.</caption><thead><tr><th><i>l</i><sub>1</sub></th><th><i>l</i><sub>2</sub></th><th><i>l</i><sub>3</sub></th></tr></thead><tbody><tr><td>6.794 · 10<sup>-5</sup></td><td>-0.01271</td><td>1</td></tr></tbody></table>

We also report in Fig.7 a comparison between the stiffness surface $K(\alpha, L)$ and the constructed point from the evaluation model with respect to the $\alpha$ and $L$ values.

The model provides a good approximation of the Evaluation Model and, as a consequence, of the experimental data, as it is shown in Fig.8. Moreover, the new model of the total force
<table><caption>Table 4: Parameters of Control Model stiffness.</caption><thead><tr><th>i</th><th>k<sub>i1</sub></th><th>k<sub>i2</sub></th><th>k<sub>i3</sub></th><th>k<sub>i4</sub></th><th>k<sub>i5</sub></th><th>k<sub>i6</sub></th><th>k<sub>i7</sub></th></tr></thead><tbody><tr><td>1</td><td>-0.00041</td><td>0.0025</td><td>-0.0061</td><td>0.00763</td><td>-0.0052</td><td>0.0018</td><td>-0.00026</td></tr><tr><td>2</td><td>0.0586</td><td>-0.352</td><td>0.8444</td><td>-1.0277</td><td>0.6697</td><td>-0.2257</td><td>0.0319</td></tr><tr><td>3</td><td>-1.67504</td><td>9.829</td><td>-22.692</td><td>25.896</td><td>-15.268</td><td>4.5087</td><td>-0.5351</td></tr><tr><td>4</td><td>-40.977</td><td>244.666</td><td>-569.9335</td><td>664.1912</td><td>-405.645</td><td>120.853</td><td>-13.045</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>Figure 7: Fitting of the stiffness K(α, L) surface with respect to α and L</figcaption></figure>

gives a good approximation of the Evaluation one in the range [0.4, 1.5]L₀, see Fig.9. This model still provides a good level of approximation of the available data. Indeed, the residual vector norm is quite small |r|₂ = 0.18, especially for L < 1.5. Since higher length values are usually not admissible in standard motions due to the mechanical limitation of biological joints, one can exploit the proposed control model in scenarios involving such standard motions.

# 5. Sensors modelling

In this section, we talk about the dynamic modeling of the organs available as sensors. They measure the muscle length and velocity, along with the applied force at the tendon level. Respectively, these organs are spindle II and Ia, and Ib (also called Golgi Tendon Organs). For a model of the last ones, we refer the reader to the works [42] and [23].

## 5.1. Spindles Ia and II

Ia- and II- spindles' neuronal activities are correlated to the length and rate of change of the muscle fibers' contractions. Along with each of the Spindles, we find active elements, i.e. the γ-motoneurons whose excitation frequency modifies the spindles' characteristics. In particular, such modulation can be interpreted as an increase in sensing precision obtained by the Spindles, despite the noise amplification drawback. We can generally distinguish between static, γ, and dynamic, γd, γ-motoneurons. In the modeling, we only consider the effects of

<figure><img src="image_10.png" /><figcaption>Figure 8: Data comparison with Spring Model. r indicates the residuals vector.</figcaption></figure>

the γd-motoneurons on the Spindles spiking frequency due to the lack of available data.

It is worth noticing in Fig.10b that in Spindle-Ia during negative velocities transient the output of the sensor saturates at zero while it is not clear what happens for Spindle-II due to the end of γs-excitation during the same part of the trajectory (negative velocities transients).

To construct the model of the sensors, we first find the static behavior of the sensors' outputs, we make an analysis on the slopes as similarly done in [24].

We now first show the hybrid dynamical model found for the Spindle-Ia and then the one of the Spindle-II.

### 5.1.1. Spindle-Ia

We considered an Hybrid dynamical model and is given by

$$ 
\begin{aligned}
C_{Ia} : \quad & \begin{cases} \dot{x}_{Ia} = p_{Ia}(K_{L_{Ia}}L - x_{Ia}) \\ \dot{p}_{Ia} = 0 \\ \dot{K}_{L_{Ia}} = 0 \\ y_{Ia} = \text{sat}_0(x_{Ia} + K_{\gamma_{Ia}}\gamma_{Ia} + K_{L_{Ia}}\dot{L}) \end{cases} \\ 
D_{Ia} : \quad & \begin{cases} (\forall X_{Ia} \in D_{Ia1}) \quad K_{L_{Ia}}^{+} = 0, \quad p_{Ia}^{+} = \bar{p}_{Ia} \\ (\forall X_{Ia} \in D_{Ia2}) \quad K_{L_{Ia}}^{+} = \bar{K}_{L_{Ia}}, \quad p_{Ia}^{+} = \bar{p}_{Ia} \\ (\forall X_{Ia} \in D_{Ia3}) \quad K_{L_{Ia}}^{+} = 0, \quad p_{Ia}^{+} = \underline{p}_{Ia} \\ (\forall X_{Ia} \in D_{Ia}) \quad x_{Ia}^{+} = x_{Ia}, \quad y_{Ia}^{+} = y_{Ia} \end{cases}
\end{aligned} 
\quad (13) 
$$

where xIa is the filter state, yIa is its output, the function sat0 is the function that saturates at zero negative values, pIa is the filter pole, KLIa, KVla and KγIa are the static gains for the muscle
<figure><img src="image_11.png" /><figcaption>Figure 9: Comparison between the Evaluation and the Control Model for different excitation frequencies, i.e. $\alpha \in 0, 30, 60, 120$ pps.</figcaption></figure>

length, velocity and the excitation state of the $\gamma_d$-motoneuron; in this case, $L, V$ and $\gamma_a$ are the muscle length, velocity and the excitation state of $\gamma_d$-motoneuron and are considered as system inputs. Then $X_{Ia}$ is the vector collecting the state, inputs and output of the system. The flow set $C_{Ia} = \mathbb{R}^3\setminus D_{Ia}$ while the jump set is $D_{Ia} = D_{Ia1} \cup D_{Ia2} \cup D_{Ia3}$, where

$$ 
\begin{aligned}
D_{Ia1} &= \{X_{Ia} : \dot{L} < 0 \& (K_{L_{Ia}} \neq 0 || p_{Ia} \neq \bar{p}_{Ia})\}, \\
D_{Ia2} &= \{X_{Ia} : \dot{L} > 0 \& (K_{L_{Ia}} \neq \bar{K}_{L_{Ia}} || p_{Ia} \neq \bar{p}_{Ia})\}, \\
D_{Ia3} &= \{X_{Ia} : \dot{L} = 0 \& (K_{L_{Ia}} \neq 0 || p_{Ia} \neq \bar{p}_{Ia})\}.
\end{aligned}
$$

In Fig.10 we compare the output of the dynamical system (13) with the available experimental data. The identified parameters are all gathered in Tab.5.

<table><caption>Table 5: Parameters of spindle Ia model.</caption><thead><tr><th>p̄<sub>Ia</sub></th><th>p<sub>Ia</sub></th><th>K̄<sub>L-Ia</sub></th><th>K<sub>L-Ia</sub></th><th>K<sub>Ia</sub></th></tr></thead><tbody><tr><td>85</td><td>3.858</td><td>7.15</td><td>5.15</td><td>0.3514</td></tr></tbody></table>

### 5.1.2. Spindle-II

Also for this sensing system a Hybrid dynamical model and is given, similarly, by

$$ 
\begin{aligned}
C_{II} : \quad & \left\{ \begin{array}{l} \dot{x}_{II} = p_{II}(\xi_{L_{II}} K_{L_{II}} L - x_{II}) \\ \dot{p}_{II} = 0 \end{array} \right. \\
& \quad \left\{ \begin{array}{l} \dot{\xi}_{L_{II}} = 0 \\ y_{II} = \text{sat}_0 \left[ x_{II} + K_{\gamma_{II}} \gamma_{II} + K_{\dot{L}_{II}} \dot{L} + (1 - \xi_{L_{II}}) K_{L_{II}} (L - \bar{L}) \right] \end{array} \right. \\
D_{II} : \quad & \left\{ \begin{array}{l} (\forall X_{II} \in D_{II1}) \quad \xi_{L_{II}}^{+} = 0, \quad p_{II}^{+} = \bar{p}_{II} \\ (\forall X_{II} \in D_{II2}) \quad \xi_{L_{II}}^{+} = 1, \quad p_{II}^{+} = \bar{p}_{II} \\ (\forall X_{II} \in D_{II3}) \quad \xi_{L_{II}}^{+} = 0, \quad p_{II}^{+} = \bar{p}_{II} \\ (\forall X_{II} \in D_{II}) \quad x_{II}^{+} = x_{II}, \quad y_{II}^{+} = y_{II} \end{array} \right.
\end{aligned}
\quad (14) $$

<figure><img src="image_12.png" /><figcaption>(a) Spindle Ia simulation with hybrid dynamical system</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(b) Experimental Data of Spindle Ia from [21][Ch.19]</figcaption></figure>

<caption>Figure 10: Spindle Ia Model output with experimental data in [21][Ch.19] comparison</caption>

where $x_{II}$ is the filter state, $y_{II}$ is system output, the function $\text{sat}_0$ is the function that saturates at zero any negative argument, $p_{II}$ is the filter pole, $K_{L_{II}}$, $K_{V_{II}}$ and $K_{\gamma_{II}}$ are the static gains for the muscle length, velocity and the excitation state of $\gamma_d$-motoneuron; in this case, $L, V$ and $\gamma_{II}$ are the muscle length, velocity and the excitation state of $\gamma_d$-motoneuron and are considered as system inputs. Then $X_{II}$ is the vector collecting the state, inputs, and output of the system. The flow set $C_{II} = \mathbb{R}^3\setminus D_{II}$ while the jump set is $D_{II} = D_{II1} \cup D_{II2} \cup D_{II3}$, where

$$ 
\begin{aligned}
D_{II1} &= \{X_{II} : \dot{L} < 0 \& (\xi_{L_{II}} \neq 0 || p_{II} \neq \bar{p}_{II})\}, \\
D_{II2} &= \{X_{II} : \dot{L} > 0 \& (\xi_{L_{II}} \neq 1 || p_{II} \neq \bar{p}_{II})\}, \\
D_{II3} &= \{X_{II} : \dot{L} = 0 \& (\xi_{L_{II}} \neq 0 || p_{II} \neq \bar{p}_{II})\}.
\end{aligned}
$$

In Fig.11 we make a comparison between the output of the dynamical system (14) with the experimental available. The identified parameters are all gathered in Tab.5.

<table><caption>Table 6: Parameters of spindle II model.</caption><thead><tr><th>p̄<sub>II</sub></th><th>p<sub>II</sub></th><th>K̄<sub>L-II</sub></th><th>K<sub>L-II</sub></th><th>K<sub>I-II</sub></th><th>K<sub>I-II</sub></th><th>L̄</th></tr></thead><tbody><tr><td>85</td><td>2.858</td><td>4.334</td><td>0.9184</td><td>0.343</td><td>2.5</td><td></td></tr></tbody></table>
<figure><img src="image_14.png" /><figcaption>Figure 11: Spindle II Model output with experimental data in [21][Ch.19] comparison</figcaption></figure>

## 6. Application Example: Oculomotor System

In this section, we provide an application of the proposed muscle modeling in the framework of the Oculomotor System. In particular, this application has been motivated by the work [43], [44], [45], in which a description of the oculomotor system has been provided with the assumption that the actuation is not given by muscular torque but rather given directly as an acceleration, see equation (11b) in [45] for example. This, in our opinion, oversimplifies the tackled problem implying that a more complex analysis must be pursued when considering an adaptive control description of the system.

We believe thus that a more insightful model has to take into account the skeletal muscle activation due to the firing rate of the α-neurons $\alpha_i$, $i = 1,2$, and that their exert force then generates an equivalent torque on the eyeball. Thus, according to the labels in Fig.12, we propose the following model for the horizontal motion of an eyeball

$$ \begin{aligned} \dot{\theta} &= \omega \\ J\dot{\omega} &= -\beta\omega + RK_1(\alpha_1, L_1)(L_1 - L_{01}(\alpha_1)) \\ &\quad - RK_2(\alpha_2, L_2)(L_2 - L_{02}(\alpha_2)) \end{aligned} \quad (15) $$

where $\theta$ is the angular displacement with respect to the sagittal

<figure><img src="image_15.png" /><figcaption>Figure 12: Schematic description of the oculomotor system in the horizontal plane evolution. The oculomotor system presents two antagonist muscles whose combined force produces a torque on the eyeball.</figcaption></figure>

(or longitudinal) plane, $\omega$ is its rate of change, $J$ is the eyeball inertia with respect to the vertical axis, $R$ is its radius, and $\alpha_i$, $K_1$, $L_i$, $L_{0i}$, for $i = 1,2$, are the α-motoneurons firing rates, the muscles variable stiffness, current and rest lengths, respectively. Moreover, the length of each muscle is related to the angle $\theta$, by considering the following relationship, $i = 1,2$,

$$ L_i = L_{i0} + (-1)^i R\theta $$

where $L_{i0}$ are the initial muscle length. And so the model can be written as a function of the system state variable

$$ \begin{aligned} \dot{\theta} &= \omega \\ J\dot{\omega} &= -\beta\omega + RK_1(\alpha_1, L_{10}-R\theta)(L_{10}-R\theta-L_{01}(\alpha_1)) \\ &\quad - RK_2(\alpha_2, L_{20}+R\theta)(L_{20}+R\theta-L_{02}(\alpha_2)) \end{aligned} \quad (16) $$

With the proposed model, we show that such an assumption does not hold any more without additional assumptions on the muscle force generation, e.g., a part of the brain is involved in providing an inversion map between the desired acceleration torque. And it is worth to notice that such an inverse map is not possible to get because the forward mapping from the firing rates $\alpha_1$ and $\alpha_2$ to the torque $\tau = RK_1(\alpha_1, L_1)(L_1 - L_{01}(\alpha_1)) - RK_2(\alpha_2, L_2)(L_2 - L_{02}(\alpha_2))$ is not injective but only surjective. Indeed, in particular, one can exploit the additional degree of freedom to obtain a desired joint stiffness for example.

Moreover, the oculomotor system involves also the sensory part described in sec.5 (one for each muscle). In particular, assuming that $L_{10} = L_{20}$ the scaled difference between $L_2$ and $L_1$ provides the angle $\theta$, while the difference between $\dot{L}_2$ and $\dot{L}_1$ provides its rate of change. These signals can be obtained from the response of the spindle Ia and II, of each muscle. Moreover, this information is in synergy with that coming from the visual cortex, when, for example pursuing an object, we also have an additional error signal reconstructed by the visual cortex itself. However, the information provided to the brain-implemented controller is asynchronous, in the sense that the time delay affecting the feedback of these signals is not the same, in general,
for the two types of information, i.e., that coming from the sensory system and that coming from the visual system.

The challenge is then to obtain a robust (to the uncertainties of the parameters) control strategy that reconstructs the system state from the Hybrid model of the Spindle *Ia* and *II* and guides the firing rates of the two muscles to achieve the desired behavior, such as object pursuing or gaze fixation to name two.

## 7. Conclusions and future works

In this work, we dealt with the modeling of muscle force considering the passive and active components of the isometric force. We first provided a brief overview of the model available in the literature and then show two models, i.e., an intermediate and a final one. The first one is more complex and precise than the second but it does not add much to the models available in the literature. On the other hand, it allowed us to construct a set of data that have been useful for the identification of the second model parameters. The latter (the spring model) is easier to handle, although less precise than the intermediate one, and gives a physical intuition about the muscle behavior. Indeed, it describes such muscular behavior as a controllable non-linear spring with controllable (variable) stiffness *K* and rest length *L*₀. We furthermore provide a hybrid dynamical model formulation of the spindles organs that act as muscle length and velocity sensors for the Central Nervous System. We then provide an application to the case of the oculomotor system proposing its dynamic modeling, and describing the related challenges.

Future works may concern the development of a muscle force characteristic depending also on the rate of shortening and lengthening of the muscle fibers' length, to provide a more complete model. Such a model can then be analyzed in a control context involving (biological) joints with antagonist actuation. Moreover, a system theoretic description of the physiological reflexes, such as the stretch reflex, can be investigated by exploiting the proposed model.

## References

[1] E. N. Marieb, K. Hoehn, Human anatomy & physiology, Pearson education, 2007.

[2] R. L. Lieber, Skeletal muscle structure, function, and plasticity, Lippincott Williams & Wilkins, 2002.

[3] W. R. Frontera, J. Ochala, Skeletal muscle: a brief review of structure and function, Calcified tissue international 96 (3) (2015) 183-195.

[4] H. L. Sweeney, D. W. Hammers, Muscle contraction, Cold Spring Harbor perspectives in biology 10 (2) (2018) a023200.

[5] I. Boyd, The isolated mammalian muscle spindle, Trends in neurosciences 3 (11) (1980) 258-265.

[6] A. V. Hill, The heat of shortening and the dynamic constants of muscle, Proceedings of the Royal Society of London. Series B-Biological Sciences 126 (843) (1938) 136-195.

[7] A. F. Huxley, Muscle structure and theories of contraction, Prog. Biophys. Biophys. Chem 7 (1957) 255-318.

[8] A. Gordon, A. F. Huxley, F. Julian, The variation in isometric tension with sarcomere length in vertebrate muscle fibres, The Journal of physiology 184 (1) (1966) 170-192.

[9] D. T. McRuer, R. E. Magdaleno, G. P. Moore, A neuromuscular actuation system model, IEEE Transactions on Man-Machine Systems 9 (3) (1968) 61-71.

[10] J. M. Winters, L. Stark, Analysis of fundamental human movement patterns through the use of in-depth antagonistic muscle models, IEEE transactions on biomedical engineering (10) (1985) 826-839.

[11] J. M. Winters, An improved muscle-reflex actuator for use in large-scale neuromusculoskeletal models, Annals of biomedical engineering 23 (4) (1995) 359-374.

[12] I. E. Brown, E. J. Cheng, G. E. Loeb, Measured and modeled properties of mammalian skeletal muscle. ii. the effectsof stimulus frequency on force-length and force-velocity relationships, Journal of Muscle Research & Cell Motility 20 (7) (1999) 627-643.

[13] G. Joyce, P. Rack, D. Westbury, The mechanical properties of cat soleus muscle during controlled lengthening and shortening movements, The Journal of physiology 204 (2) (1969) 461-474.

[14] P. M. Rack, D. Westbury, The effects of length and stimulus rate on tension in the isometric cat soleus muscle, The Journal of physiology 204 (2) (1969) 443-460.

[15] J. M. Winters, Hill-based muscle models: a systems engineering perspective, in: Multiple muscle systems, Springer, 1990, pp. 69-93.

[16] F. C. Van der Helm, L. A. Rozendaal, Musculoskeletal systems with intrinsic and proprioceptive feedback, in: Biomechanics and neural control of posture and movement, Springer, 2000, pp. 164-174.

[17] N. Hogan, Adaptive control of mechanical impedance by coactivation of antagonist muscles, IEEE Transactions on automatic control 29 (8) (1984) 681-690.

[18] J. McIntyre, J.-J. E. Slotine, Does the brain make waves to improve stability?, Brain research bulletin 75 (6) (2008) 717-722.

[19] J. M. Wakeling, S. S. Lee, A. S. Arnold, M. de Boef Miara, A. A. Biewener, A muscle's force depends on the recruitment patterns of its fibers, Annals of biomedical engineering 40 (8) (2012) 1708-1720.

[20] I. Boyd, P. Murphy, V. Moss, Analysis of primary and secondary afferent responses to stretch during activation of the dynamic bag 1 fibre or the static bag 2 fibre in cat muscle spindles, in: The muscle spindle, Springer, 1985, pp. 153-158.

[21] I. A. Boyd, M. H. Gladden, The muscle spindle, Springer, 1985.

[22] U. Windhorst, Muscle proprioceptive feedback and spinal networks, Brain research bulletin 73 (4-6) (2007) 155-202.

[23] M. P. Mileusnic, G. E. Loeb, Mathematical models of proprioceptors. ii. structure and function of the golgi tendon organ, Journal of neurophysiology 96 (4) (2006) 1789-1802.

[24] M. G. Maltenfort, R. Burke, Spindle model responsive to mixed fusimotor inputs and testable predictions of β feedback effects, Journal of neurophysiology 89 (5) (2003) 2797-2809.

[25] E. Schomburg, Spinal sensorimotor systems and their supraspinal control, Neuroscience research 7 (4) (1990) 265-340.

[26] P. B. Matthews, The human stretch reflex and the motor cortex, Trends in neurosciences 14 (3) (1991) 87-91.

[27] E. Jankowska, A neuronal system of movement control via muscle spindle secondaries., Progress in brain research 80 (1989) 299-303.

[28] M. M. Gassel, A critical review of evidence concerning long-loop reflexes excited by muscle afferents in man, Journal of Neurology, Neurosurgery & Psychiatry 33 (3) (1970) 358-362.

[29] R. M. Eccles, A. Lundberg, Supraspinal control of interneurones mediating spinal reflexes, The Journal of physiology 147 (3) (1959) 565-584.

[30] K. Iqbal, A. Roy, Stabilizing pid controllers for a single-link biomechanical model with position, velocity, and force feedback, J. Biomech. Eng. 126 (6) (2004) 838-843.

[31] D. A. Kistemaker, A. J. K. Van Soest, J. D. Wong, I. Kurtzer, P. L. Gribble, Control of position and movement is simplified by combined muscle spindle and golgi tendon organ feedback, Journal of neurophysiology 109 (4) (2013) 1126-1139.

[32] E. Mtui, G. Gruener, P. Dockery, Fitzgerald's Clinical Neuroanatomy and Neuroscience E-Book, Elsevier Health Sciences, 2020.

[33] J. Nolte, The human brain, Mosby/Elsevier, 1993.

[34] D. M. Wolpert, R. C. Miall, M. Kawato, Internal models in the cerebellum, Trends in cognitive sciences 2 (9) (1998) 338-347.

[35] O. Barak, Mark I. latash, neurophysiological basis of movement. human kinetics, champaign, 1998, Medicinski pregled 55 (5-6) (2002) 261-261.

[36] J. A. Kiernan, M. L. Barr, Barr's the human nervous system: an anatomical viewpoint, Lippincott Williams & Wilkins, 2009.

[37] M. A. Patestas, L. P. Gartner, A textbook of neuroanatomy, John Wiley & Sons, 2016.
[38] J. Huang, A. Isidori, L. Marconi, M. Mischiati, E. Sontag, W. Wonham, Internal models in control, biology and neuroscience, in: 2018 IEEE Conference on Decision and Control (CDC), IEEE, 2018, pp. 5370–5390.

[39] D. M. Wolpert, M. Kawato, Multiple paired forward and inverse models for motor control, Neural networks 11 (7-8) (1998) 1317–1329.

[40] M. Kawato, Internal models for motor control and trajectory planning, Current opinion in neurobiology 9 (6) (1999) 718–727.

[41] M. Mischiati, H.-T. Lin, P. Herold, E. Imler, R. Olberg, A. Leonardo, Internal models direct dragonfly interception steering, Nature 517 (7534) (2015) 333–338.

[42] M. P. Mileusnic, I. E. Brown, N. Lan, G. E. Loeb, Mathematical models of proprioceptors. i. control and transduction in the muscle spindle, Journal of neurophysiology 96 (4) (2006) 1772–1788.

[43] M. E. Broucke, Model of the oculomotor system based on adaptive internal models, IFAC-PapersOnLine 53 (2) (2020) 16430–16437.

[44] M. E. Broucke, Adaptive internal model theory of the oculomotor system and the cerebellum, IEEE Transactions on Automatic Control 66 (11) (2020) 5444–5450.

[45] E. Battle, M. E. Broucke, Adaptive internal models in the optokinetic system, in: 2021 60th IEEE Conference on Decision and Control (CDC), IEEE, 2021, pp. 641–648.