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

The muscle spindle model presented in Maltenfort and Burke (2003) calculates muscle spindle primary aﬀerent feedback depending on the muscle ﬁbre stretch and fusimotor drive. The aim of this paper is to provide an updated version of the model, which is now capable of replicating the originally published data. This is achieved by modifying the equations describing the modulation of the muscle spindle output in response to dynamic fusimotor drive.


Introduction
Muscle spindles can be found in almost all skeletal muscles and make the most important contribution to proprioception (Macefield and 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 α-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 (β )-as well as static and dynamic fusimotor (γ)-neurons (Banks, 1994;Matthews, 1962). The muscle spindle model by Maltenfort and 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 and Matthews, 1964;Hulliger et al., 1977a,b). 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 and Maier (2014) previously adapted the Maltenfort and 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 and 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 and Burke (2003) as original model and to our model as the updated model.

Model description
The Maltenfort and Burke muscle spindle model was developed for two purposes: the investigation of β -feedback loops (Burke and 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 and Thoden (1968a,b) 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 and 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 β -innervation is still under debate, the β -feedback loop is not part of the provided model.

Model 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 velocitydependent terms modulated by dynamic fusimotor drive Q d and S 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.

General model modifications
Eqn. 2A of the original paper models the muscle spindle Ia firing rate R as the sum of a passive contribution R passive and the mixed contribution from static (∆R static ) and dynamic (∆R dynamic ) fusimotor drive: Thereby, f occlusion is a function describing the non-linear summation of the contributions from static and dynamic fusimotor input.

2/8
In detail, the contribution modulated by dynamic fusimotor drive ∆R dynamic is described by a function depending on dynamic gamma drive γ 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 Note that only the descriptions of a, c and d are modified with respect to the original publication.

Implementation 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 lag (cf. Eqn. 3C, original paper): Thereby, τ is a time constant, specific for each type of fusimotor drive. In the original publication, x 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. 9) 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 lag (t = 0) = 0 is used.

MATLAB
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

CellML
The CellML simulations were run with OpenCOR version 0.6 (Garny and 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).

Reproducability 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  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).

Discussion
With this publication, a corrected version of the Maltenfort and 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 and 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 and 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 and 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.

6/8
There are details/documentation on how the source code was compiled There are details on how to run the code in the provided documentation The initial conditions are provided for each of the simulations Details for creating reported graphical results from the simulation results Source code: a declarative language is used (e.g. SBML, CellML, NeuroML) The algorithms used are defined or cited in previous articles The algorithm parameters are defined Post-processing of the results are described in sufficient detail

Executable model provided:
The model is executable without source (e.g. desktop application, compiled code, online service) There are sufficient details to repeat the required simulation experiments

The model is described mathematically in the article(s):
Equations representing the biological system There are tables or lists of parameter values There are tables or lists of initial conditions

Machine-readable tables of initial conditions
The simulation experiments using the model are described mathematically in the article: Integration algorithms used are defined Stochastic algorithms used are defined Random number generator algorithms used are defined Parameter fitting algorithms are defined The paper indicates how the algorithms yield the desired output