Spindle Model Responsive to Mixed Fusimotor Inputs: an updated version of the Maltenfort and Burke (2003) model

1Introduction¶

Muscle spindles can be found in almost all skeletal muscles and make the most important contribution to proprioception Macefield & Knellwolf, 2018. The muscle spindle is a special mechano-receptor sensing the intrafusal muscle fiber length change and providing stretch feedback to the neuromuscular system via two types of axons: primary (Ia) and secondary (II) fibers. The afferent axons form feedback loops exciting $\alpha$-motoneurons and thus, ultimately control skeletal muscle contraction Kandel et al., 2000. For optimizing motor control in various conditions, the sensitivity of muscle spindles is modulated by specialized neurons, i. e. skeletofusimotor ($\beta$)- as well as static and dynamic fusimotor ($\gamma$)-neurons Banks, 1994Matthews, 1962.

The muscle spindle model by Maltenfort & Burke (2003) predicts primary afferent feedback depending on the muscle stretch and fusimotor drive. The model was validated using data recorded from cat muscle spindles during ramp-and-hold and sinusoidal stretches Crowe & Matthews, 1964Hulliger et al., 1977Hulliger et al., 1977.

The strength of the model is its capability to predict muscle spindle primary activity over a large range of physiological conditions while demanding little computational resources.

However, in the presence of dynamic fusimotor drive, the published model predicts spindle frequencies that are considerably higher than shown in the original publication. Grandjean & Maier (2014) previously adapted the Maltenfort & Burke (2003) model to improve the predicted physiological behaviour in the presence of fusimotor drive. However, they did not publish their code and the shown results could not be reproduced with the provided equations. Thus, the aim of this manuscript is to provide an updated and validated version of the muscle spindle model originally published in Maltenfort & Burke (2003). In detail, this work presents a model capable of replicating the muscle spindle firing rates during ramp-and-hold and sinusoidal stretches while considering variable fusimotor drives as shown in Fig. 4 and Fig. 5 of the original publication. Since the published equations are based on a specific discretization scheme and hence, the time step size is a variable, the model is implemented in MATLAB. Nevertheless, to promote the use of open standards for model sharing, we additionally provide a CellML implementation of the updated model.

Note that we will refer to the model based on the equations published in Maltenfort & Burke (2003) as original model and to our model as the updated model.

2Model description¶

The Maltenfort and Burke muscle spindle model was developed for two purposes: the investigation of $\beta$-feedback loops Burke & Tsairis, 1977 and for serving as computationally efficient model in large-scale simulations of the neuromuscular system. The latter is ensured by employing only algebraic equations. The model has a block-wise structure (cf. Fig. 1B in the original manuscript), calculating a baseline muscle spindle frequency for passive length changes as well as additional contributions modulated by static and dynamic fusimotor drives. Following the findings from Andersson et al. (1968), Lennerstrand & Thoden (1968)Lennerstrand & Thoden (1968) and Lennerstrand (1968) and using a phenomenological approach, each contribution is calculated as the sum of four components: a pure velocity and a pure position sensitivity, a mixed velocity and position sensitivity and a baseline firing rate at the initial length. Thereby, the velocity is the filtered derivative of the the position input to achieve smooth transitions from shortening to lengthening cf. Maltenfort & Burke, 2003, Eqn. 3. The static and dynamic contributions are non-linearly mixed by an occlusion function, partially suppressing the smaller contribution as it was observed in experimental studies Schäfer, 1974. Finally, the passive contribution and the mixed fusimotor contribution are added up yielding the firing rate of the muscle spindle primary axon.

Since the role of $\beta$-innervation is still under debate, the $\beta$-feedback loop is not part of the provided model.

3Model modifications¶

Comparisons with the original code, which was kindly provided by the author, revealed deviations with respect to the published equations. This exclusively concerns the descriptions of the velocity-dependent terms modulated by dynamic fusimotor drive $Q_{\mathrm d}$ and $S_\mathrm{v,d}$ from Table 1 in the original manuscript. For better comprehensibility of the quantities and their context, we will display further equations. The parameter values are all adopted from the original publication.

Note that the modifications introduced in Section 3.1 apply to both provided implementations, i. e. MATLAB and CellML, while the modifications introduced in Section 3.2 only apply to the CellML implementation.

3.1General model modifications¶

Eqn. 2A of the original paper models the muscle spindle Ia firing rate $R$ as the sum of a passive contribution $R_{\mathrm{passive}}$ and the mixed contribution from static ($\Delta R_{\mathrm{static}}$) and dynamic ($\Delta R_{\mathrm{dynamic}}$) fusimotor drive:

Thereby, $f_{\mathrm{occlusion}}$ is a function describing the non-linear summation of the contributions from static and dynamic fusimotor input.

In detail, the contribution modulated by dynamic fusimotor drive $\Delta R_{\mathrm{dynamic}}$ is described by a function depending on dynamic gamma drive $\gamma_\mathrm{d}$, displacement $x$ and velocity $v$ (cf. Eqn. 2C, original publication):

The individual terms in Eqn. 2 are specified in Table 1 as well as Eqn. 5B and 6B of the original publication. The modifications made for this publication concern the calculation of $S_\mathrm{v,d}(\gamma_\mathrm{d},v)$ and $Q_\mathrm{d}(\gamma_\mathrm{d},v)$, i.e.

Note that only the descriptions of $a$, $c$ and $d$ are modified with respect to the original publication.

3.2Implementation in CellML¶

Implementing the model in CellML requires further modifications concerning the velocity filter, which smoothes transitions between ramp stretch and shortening. In detail, the velocity $v$ is obtained from the temporally filtered position input $x_\mathrm{lag}$ (cf. Eqn. 3C, original paper):

Thereby, $\tau$ is a time constant, specific for each type of fusimotor drive. In the original publication, $x_\mathrm{lag}$ is calculated by a digital filter function (cf. Eqn. 3B, original paper) using the chosen time step size as a variable. This can not be implemented in CellML in a straight forward way. Thus, the CellML code numerically integrates the differential equation (Eqn. 4) to provide the actual velocity $v(t)$. Note that in the MATLAB code, the same formulation as in the original paper is used. In both codes, the initial condition $x_\mathrm{lag}(t=0) = 0$ is used.

4Computational simulation¶

4.1MATLAB¶

The simulations were run with MATLAB (The MathWorks, Inc., Natick, Massachusetts, United States) Version R2018a (9.4.0.813654) using a time step size of 1 ms. Scripts are provided to run the simulations and plot the results shown in Fig. 1, Fig. 2 and Fig. 3:

• Fig01.m

• Fig02.m

• Fig03.m

These scripts call the function muscle_spindle_maltenfort_and_burke_2003.m, which contains the muscle spindle model equations. Note that the figure scripts contain a boolean variable plot_paper_data, which determines whether the digitized data from the original publication (cf. Section 5) is plotted for comparison. If the variable is set true, the files

• Maltenfort_and_Burke_2003_Fig4A_digitized.csv

• Maltenfort_and_Burke_2003_Fig4B_digitized.csv

• Maltenfort_and_Burke_2003_Fig5A_digitized.csv

• Maltenfort_and_Burke_2003_Fig5B_digitized.csv

• Maltenfort_and_Burke_2003_Fig5C_digitized.csv

• Maltenfort_and_Burke_2003_Fig5D_digitized.csv

must be downloaded to the same folder as the respective figure script, such that the data can be loaded for plotting. Further, to add the results from CellML simulations to Fig. 3, the CellML simulations results need to be stored to a file called Fig03.csv and the respective boolean variable in the MATLAB script Fig03.m (plot_cellml_data) must be set true.

4.2CellML¶

The CellML simulations were run with OpenCOR version 0.6 Garny & Hunter, 2015, Euler forward solver and a time step size of 0.1 ms. The CellML simulation results are obtained by running Fig03.sedml. Thereby, the stretch scenarios shown in Figure 3A and B are applied consecutively. Note that in Figure 3B the second of two sinusoidal stretch cycles is plotted to exclude effects caused by the initial displacement. The comparisons to the respective MATLAB results, as also shown in Figure 3, are created by storing the CellML simulation results to a file called Fig03.csv and running Fig03.m (cf. Section 4.1).

5Reproducibility goals¶

The aim of this paper is to provide a model which is able to reproduce the results as shown in Fig. 4A and Fig.4B as well as Fig. 5A-D of the original publication. Since the data from the original publication were not available, Engauge Digitizer (Version 10.4) was used to obtain the data points shown herein.

In Fig. 1 spindle responses to ramp-and-hold stretches of two different ramp velocities are shown. The prediction of the simulations coincide very well with the digitized data from the original publication. The same applies for simulation results in response to sinusoidal stretches as shown in Fig. 2. Here, the second of two stretch cycles is plotted to exclude effects from the initial displacement at the beginning of the first stretch cycle. Note that small deviations from the digitized data have their origin in unavoidable inaccuracies due to figure resolution and choice of line styles.

The CellML code was run exemplary for two scenarios to show the consistency with the MATLAB model. As it can be seen in Fig. 3A, the CellML prediction is slightly higher during the ramp phase of the ramp-and-hold scenario. This is due to the necessary modification to the implementation of the velocity filter (cf. Section 3) and the difference is dependent on the chosen solver and time step size. However, the deviations are small and even less pronounced for the sinusoidal stretch with applied fusimotor drive (cf. Fig. 3B).

6Discussion¶

With this publication, a corrected version of the Maltenfort & Burke (2003) muscle spindle model is made available to facilitate the integration of muscle spindle models in future models of the neuromuscular system. The overestimation of the spindle firing rate in response to dynamic fusimotor drive, as also addressed by Grandjean & Maier (2014), is corrected in the provided model. It was shown that the provided codes reproduce the results published in the original paper. As shown in Maltenfort & Burke (2003), these results match a variety of experimentally observed characteristics of muscle spindle firing. For an elaborate discussion of the model, the reader is referred to the original publication Maltenfort & Burke, 2003.

In absence of secondary afferent firing, the model is well suited for large-scale bio-physically motivated simulations due to its low computational cost cf. e. g. Röhrle et al., 2019.

For example, combining the spindle model with a continuum-mechanical multi-muscle model e.g. Röhrle et al., 2017 has the potential to provide new means of investigating muscle spindle behaviour in response to heterogeneous muscle deformations as well as sensory interactions across joints.