An Approach for Simulation of the Muscle Force Modeling It by Summation of Motor Unit Contraction Forces

Muscle force is due to the cumulative effect of repetitively contracting motor units (MUs). To simulate the contribution of each MU to whole muscle force, an approach implemented in a novel computer program is proposed. The individual contraction of an MU (the twitch) is modeled by a 6-parameter analytical function previously proposed; the force of one MU is a sum of its contractions due to an applied stimulation pattern, and the muscle force is the sum of the active MUs. The number of MUs, the number of slow, fast-fatigue-resistant, and fast-fatigable MUs, and their six parameters as well as a file with stimulation patterns for each MU are inputs for the developed software. Different muscles and different firing patterns can be simulated changing the input data. The functionality of the program is illustrated with a model consisting of 30 MUs of rat medial gastrocnemius muscle. The twitches of these MUs were experimentally measured and modeled. The forces of the MUs and of the whole muscle were simulated using different stimulation patterns that included different regular, irregular, synchronous, and asynchronous firing patterns of MUs. The size principle of MUs for recruitment and derecruitment was also demonstrated using different stimulation paradigms

A General Mathematical Algorithm for Predicting the Course of Unfused Tetanic Contractions of Motor Units in Rat Muscle

<div><p>An unfused tetanus of a motor unit (MU) evoked by a train of pulses at variable interpulse intervals is the sum of non-equal twitch-like responses to these stimuli. A tool for a precise prediction of these successive contractions for MUs of different physiological types with different contractile properties is crucial for modeling the whole muscle behavior during various types of activity. The aim of this paper is to develop such a general mathematical algorithm for the MUs of the medial gastrocnemius muscle of rats. For this purpose, tetanic curves recorded for 30 MUs (10 slow, 10 fast fatigue-resistant and 10 fast fatigable) were mathematically decomposed into twitch-like contractions. Each contraction was modeled by the previously proposed 6-parameter analytical function, and the analysis of these six parameters allowed us to develop a prediction algorithm based on the following input data: parameters of the initial twitch, the maximum force of a MU and the series of pulses. Linear relationship was found between the normalized amplitudes of the successive contractions and the remainder between the actual force levels at which the contraction started and the maximum tetanic force. The normalization was made according to the amplitude of the first decomposed twitch. However, the respective approximation lines had different specific angles with respect to the ordinate. These angles had different and non-overlapping ranges for slow and fast MUs. A sensitivity analysis concerning this slope was performed and the dependence between the angles and the maximal fused tetanic force normalized to the amplitude of the first contraction was approximated by a power function. The normalized MU contraction and half-relaxation times were approximated by linear functions depending on the normalized actual force levels at which each contraction starts. The normalization was made according to the contraction time of the first contraction. The actual force levels were calculated initially from the recorded tetanic curves and subsequently from the modeled curves obtained from the summation of all models of the preceding contractions (the so called “full prediction”). The preciseness of the prediction was verified by two coefficients estimating the error between the modeled and the experimentally recorded curves. The proposed approach was tested for 30 MUs from the database and for three additional MUs, not included in the initial set. It was concluded that this general algorithm can be successfully used for modeling of a unfused tetanus course of a single MU of fast and slow type.</p></div

Approximation of the relationships between the angles <i>α</i>^<i>(j)</i> and the parameter <i>F_mftf^(j)</i> (reflecting the maximum force that the respective MU can develop in the fused tetanus) normalized to the amplitude of the first decomposed contraction.

<p>The two parameters used to calculate data presented on the abscissa are illustrated in a frame below the axis on recordings of a fragment of an unfused tetanus and the fused tetanus (red lines indicate amplitudes of the first twitch and the maximum tetanus force). The data for the angles are given in the fourth column of <a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0162385#pone.0162385.t001" target="_blank">Table 1</a>. S MU—blue asterisks; FR MU—red asterisks; FF—green asterisks. The black dashed curves are different approximations: ‘o’—with a linear model: <i>y = ax+b</i>, <i>a</i> = -2.994, <i>b</i> = 91.96; ‘◊‘—with an exponential model from 1^st type: <i>y = ae</i>^<i>bx</i>, <i>a</i> = 95.8, <i>b</i> = -0.04602; ‘✻‘—with an exponential model from 2^nd type: <i>y = ae</i>^<i>bx</i><i>+ce</i>^<i>dx</i>, <i>a</i> = 66.53, <i>b</i> = -0.2896, <i>c</i> = 54.6, <i>d</i> = 0.001297; ‘□’—with a power model (this model is chosen for further modeling and is marked with the bold dashed line): <i>y = ax</i>^<i>b</i>, <i>a</i> = 108.8, <i>b</i> = -0.2603. Here, <i>y = α1</i>^<i>(j</i>) and <i>x = F</i>_<i>mftf</i>^<i>(j)</i><i>/F</i>_<i>max</i>^<i>(j)</i><i>(1)</i>.</p

Dependencies between two normalized parameters: <i>F_max^(j)(i)/F_max^(j)(1)</i> and <i>F_res^(j)(i)/F_max^(j)(1)</i>, for all 30 MUs.

<p>The two parameters used to calculate data presented on the ordinate are illustrated in a frame left to the axis on an example of a train of decomposed twitches (red lines indicate amplitudes of the first and the <i>i</i>-th twitch). Additionally, the parameters used to calculate data presented on the abscissa are illustrated in a frame below the axis on a fragment of the unfused tetanus and the fused tetanus recordings (red lines indicate amplitudes of the first twitch and the residual force for the response to the <i>i</i>-th stimulus). The symbols on the main chart marked in blue present the data for slow MUs, in red—the data for FR MUs, and in green—the data for FF MUs. The data for each MU was approximated by straight lines in respective colors: blue for S MUs, red for FR MUs, and green for FF MUs. The angles <i>α</i>^<i>(j)</i> (<i>j</i> = 1, 2,…, 30) for each MU were calculated between these lines and the ordinate.</p

Comparison between experimental tetani and the force curves, obtained through the new approximation approach.

<p>Full prediction of three additional tetanic curves, obtained by stimulation with three new patterns, applied to three MUs not included in the input database. Red color—the recorded curve; blue color—the predicted force curve. A. A slow MU stimulated with the mean frequency of 20 Hz (S11 in <a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0162385#pone.0162385.t001" target="_blank">Table 1</a>), the stimulation pattern with interpulse intervals IPI4 is used; B. A FR MU stimulated with the mean frequency of 50 Hz (FR11 in <a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0162385#pone.0162385.t001" target="_blank">Table 1</a>), the stimulation pattern with interpulse intervals IPI5 is used; C. A FF MU stimulated with the mean frequency of 33.3 Hz (FF11 in <a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0162385#pone.0162385.t001" target="_blank">Table 1</a>), the stimulation pattern with interpulse interval IPI6 is used. The values of the angles <i>α1</i>^<i>(i)</i> are shown in <a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0162385#pone.0162385.t001" target="_blank">Table 1</a>. IPIs are given in <a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0162385#pone.0162385.s001" target="_blank">S1 File</a> Table B. Note that the time and the force scales are different for the three MUs.</p

Illustration of the calculation of the parameter <i>F^pr_tetmin^(j)</i>(<i>i</i>) for the 23^th MU after adding the models of all preceding contractions.

<p>The first five models (plotted by using black dashed lines) are summed and the accumulated force is plotted as a solid black line; the value of the force <i>F</i>^<i>pr</i>_{<i>tetmin</i>}^<i>(23</i>) (6) is computed at the moment when the 6^th pulse comes and the next, 6^th contraction, is calculated using this value. The model of this 6^th contraction is red dashed line. The solid black line after the appearance of the 6-th pulse shows the addition of the force evoked by the 6^th pulse to the previous five contractions.</p

Description of the parameters used for the <i>j</i>-th MU recordings.

<p>A. The model of the <i>i</i>-th decomposed contraction within the unfused tetanus is shown by a dashed line, the black solid line is a piece of the force obtained by subtraction of all previous (<i>i</i>-1) contraction models from the experimental tetanic force. The parameters of the <i>i</i>-th twitch-like contraction are: <i>F</i>_<i>max</i>^<i>(j)</i><i>(i)</i>–the maximum twitch force; <i>T</i>_<i>lead</i>^<i>(j)</i> <i>(i)</i>–the lead time, the time between the <i>i</i>-th stimulus (its time position is indicated by vertical arrow) and the start of the current <i>i</i>-th contraction; <i>T</i>_<i>hc</i>^<i>(j)</i>–the half-contraction time, the time from the start of the contraction to the moment when the twitch force reaches a half of its maximal value; <i>T</i>_<i>c</i>^<i>(j)</i>–the contraction time, the time from the start of the contraction to the moment when the twitch amplitude reaches it maximal value <i>F</i>_<i>max</i>^<i>(j)</i><i>(i); T</i>_<i>hr</i>^<i>(j)</i><i>(i)</i>–the half-relaxation time, the time between the start of the contraction to the moment when during the relaxation, the twitch force decreases to <i>F</i>_<i>max</i>^<i>(j)</i><i>(i)/</i>2; <i>T</i>_<i>tw</i>^<i>(j)</i><i>(i)</i>—the duration of the current contraction, from the time between the moment when the contraction starts and the moment when the force decreases to 0.01% of <i>F</i>_<i>max</i>^<i>(j)</i><i>(i)</i>. The equation describing this bell-shape 6-parameters curve are given in [<a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0162385#pone.0162385.ref009" target="_blank">9</a>]; B. Parameters measured for tetanic contractions presented on a part of the unfused tetanic curve (left) and the maximum fused tetanus (right). <i>F</i>_<i>mftf</i>^<i>(j)</i>—the maximal force that a MU develops during stimulation at 150 Hz stimulation frequency (in the fused tetanus). <i>F</i>_{<i>tetmin</i>}^<i>(j)</i><i>(i)</i>—the force level at which the <i>i</i>-th contraction starts; <i>F</i>_<i>res</i>^<i>(j)</i><i>(i) = F</i>_<i>mftf</i>^<i>(j)</i><i>-F</i>_{<i>tetmin</i>}^<i>(j)</i><i>(i)</i>—the residual force.</p

Calculated angles and coefficients reflecting the similarity between the experimental curves and the predicted curves for all 33 MUs.

<p><i>FitCo1</i> and <i>AreaCo1</i> are these coefficients when the modeled curve is obtained as a sum of equal to the model of the first contraction twitches, according to the respective stimulation pattern. <i>α1</i>^<i>(j)</i> is the angle calculated by using the Eq (<a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0162385#pone.0162385.e003" target="_blank">1</a>). <i>FitCo2</i> and <i>AreaCo2</i> are the coefficients for the experimental and modeled curve (the calculated experimental values of <i>F</i>_{<i>tetmin</i>} ^<i>(j)</i><i>(i)</i> are used as input parameters for the prediction). <i>FitCo3</i> and <i>AreaCo3</i> are the coefficients for the experimental and modeled curve using the same angles <i>α1</i>^<i>(j)</i> but for full prediction algorithm (i.e. values of F_tetmin ^(j)(i) are also predicted consecutively). <i>α2</i>^<i>(j)</i> is the angle obtained by a sensitivity analysis so that the modeled curve is the most similar to the experimental one, and <i>FitCo4</i> and <i>AreaCo4</i> are the respective coefficients. <i>α3</i>^<i>(j)</i> is the improved angle using the Eq (<a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0162385#pone.0162385.e013" target="_blank">11</a>) (see <a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0162385#pone.0162385.g008" target="_blank">Fig 8</a>), and <i>FitCo5</i> and <i>AreaCo5</i> are the respective coefficients.</p