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

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.


Introduction
The mechanical output of a skeletal muscle is the sum of forces generated by its active motor units (MUs), and the forces of the individual MUs depend on a firing pattern of the motoneurons. During voluntary activity, motoneurons generate trains of action potentials at variable time intervals [1][2][3][4][5], and each of these action potentials evokes a twitch-like force response [6][7][8]. However, when MU tetanic contractions evoked by trains of stimuli at variable time intervals (which cause considerable force fluctuations) are decomposed into twitch-shape responses to successive stimuli, the amplitude and the time parameters of these responses appear to be highly variable, and considerable differences in the tetanic force development are observed between the three types of MUs [9,10]. For slow (S) MUs, individual twitch-like responses have force amplitudes up to seven times higher and contraction and relaxation phases over three times longer in comparison to the single twitch [11]. For fast resistant-to-fatigue (FR) MUs, variability of the twitch parameters is in general much lower (force amplitudes up to three times higher than the single twitch), whereas for some fast fatigable (FF) MUs the decomposed twitch-like responses frequently have even lower amplitudes than the amplitude of the single twitch. The parameters of the twitch-like responses to individual stimuli appear to depend mainly on the force level reached by a MU when the next stimulus is delivered. This dependence is specific for different MU types [10].
Muscle models composed of MUs [12,13,14] give deeper and new insight into the processes of generation and control of the muscle force. In each model, there are many simplifications and assumptions that can lead to errors influencing the muscle model correctness. Summation of equal twitches [15] is the easiest, the most applicable and logical method for modeling tetanic contractions. However, the error made by using this method seems essential, especially for slow MUs [10,15]. The principal questions are: how do the successive contractions change in comparison to the single twitch? Which parameters determine these changes? How can they be modeled? In the study by Fuglevand et al. (1993) [12] a gain was used for changing the amplitude of successive contractions and this parameter depended on interpulse intervals (IPIs). However, this gain was not specific for the MU type. Moreover, our previous studies have demonstrated that contraction as well as relaxation time parameters of the successive contractions are also variable [10], and that the actual force level is a better predictor of the tetanic force increase than the IPIs [16]. A suitable, physiologicallybased approach for prediction of these phenomena, differentiating the three MU types, is still missing.
The aim of the present paper was to find a general rule for determining the force development for various types of MUs and to develop a new, physiologically-based, mathematical approach for predicting the parameters of twitch-like responses to successive stimuli with different IPIs, delivered in random patterns, which evoke moderately fused tetanic contractions. The derived equations were based on decomposition of tetanic curves of 30 MUs (10 from each physiological type). In order to validate and verify accuracy of the developed mathematical approach we accomplished the following: (1) predicted the profiles of unfused tetanic contractions for the 30 MUs from the input database used to develop the model; (2) compared the modeled force profiles to the recordings from physiological experiments and to contractions obtained by summation of equal twitches; (3) applied the proposed general mathematical algorithm for prediction of tetanic contractions of three additional MUs outside the database, and not processed during the decomposition and modeling, furthermore evoked at three different patterns of stimuli.

Experimental data
The input data for this modeling study were taken from the results of the decomposition of the recorded curves of unfused tetanic contractions for 30 MUs of the rat medial gastrocnemius muscle- Fig 1. We decomposed 10 MUs of each type: slow (S), fast resistant to fatigue (FR) and fast fatigable (FF). The electrophysiological experiments for the recorded tetanic force curves were presented in our previous paper [11]. The experimental procedures (it was approved by the Local Bioethics Committee in Poznan-permission No 2/2015) of in vivo recordings from functionally isolated MUs were there explained in detail. Most importantly, for each MU the following parameters were recorded: the single twitch (Fig 1A), the maximum tetanus at 150 Hz constant rate stimulation (the fused tetanus- Fig 1C) and several unfused tetanic contractions ( Fig 1B) at stimulation frequencies evoking tetanus with force levels ranging from 30 to 70% of the maximum tetanus force, and with fusion indices in the range of 0.40 to 0.95. For definitions of these parameters). B. The initial part of an unfused tetanic contraction recording (blue) and models of the successive contractions (black) obtained using a mathematical decomposition of the recording into force responses to individual stimuli. Green dots indicate force levels at the beginnings of the twitchshape responses to successive stimuli. The red dots indicate the stimuli. C. A recording of a fused tetanic contraction. The red vertical line indicates the maximum tetanus force. D. The amplitudes of the decomposed twitches as a function of the force levels developed by a motor unit at the time moments when the contractions begin. The main aim of the modelling is to predict the course of an unfused tetanic contraction (blue recording) developed during stimulation at variable interpulse intervals using as input data the parameters indicated in red, i.e. the twitch time and force parameters, the maximum tetanus force and the applied pattern of stimuli. The previously observed correlation between the force of the decomposed twitches and the present MUs' force level gave us a start point for the modelling.
doi:10.1371/journal.pone.0162385.g001 each frequency the unfused tetanus was evoked by two stimulation patterns composed of 41 pulses. Initially, stimulation at a constant frequency with equal IPIs was applied, and afterwards, stimulation at the same mean frequency, but with non-equal IPIs and within a range of 50-150% of the mean IPI was used [13,16]. The used MU firing frequencies were: from 10 Hz to 16

Decomposition of the unfused tetanic curve
Despite the fact that several patterns, which evoked contractions with considerably variable fusion degree, were recorded using different stimulation frequencies, for each of the chosen 30 MUs only one middle-fused tetanic curve with non-equal IPIs was taken for further decomposition and statistics [11]. We used 13 different stimulation patterns of the following frequencies: one of 50 Hz, three of 40 Hz, one of 33.3 Hz, one of 25 Hz, one of 20 Hz, two of 16.6 Hz, two of 14.3 Hz, one of 12.5 Hz and one of 10 Hz (6 patterns are given in S1 File Table B). The experimental curves and the respective IPIs were recorded as text files. These 30 tetanic curves were decomposed using the method described previously [17]. The decomposed 41 successive twitch-like contractions ( Fig 1B) were modeled using the 6-parameter analytical function proposed and verified in [9]. The data generalizing all calculated parameters, i.e. the 6 parameters of all decomposed contractions and the maximal force of the fused tetanus, were stored for further analysis.
To verify the accuracy of the decomposition process, a tetanic curve was modeled for each MU, calculating the sum of the models of all decomposed twitch-like contractions according to time distribution of successive pulses in the applied stimulation pattern. The similarity between the recorded and the modeled curves was estimated by using two coefficients proposed and used in [15]. The first one was the fit coefficient, calculated using the formula where N is the number of samples and Δ i is the i-th difference between the two curves, normalized to their common maximal value. When the recorded and the modeled curves match perfectly, then FitCo = 100%; the lower the FitCo, the bigger the difference between them. The other coefficient, AreaCo, compares the area under the curves. If area exp and area mod are the areas under the recorded and the modeled curves, respectively, AreaCo ¼ area exp area mod . If the value of AreaCo is 1, the modeled and the recorded forces are similar; if it is higher than 1, the recorded force output of this MU is bigger. If it is less than 1, the model predicts higher total force output during the respective MU activity.
To estimate the error that could be made using a simple algebraic summation of equal twitches (i.e. equal to the model of the first decomposed contraction) using the same stimulation pattern, the same two coefficients were calculated for the experimental and the so modeled curves.

Approximation of the data obtained from the decomposition of the tetanic curves
The data from the decomposition of the 30 tetanic curves of the chosen MUs were further used in the approximation procedure. Each decomposed contraction was modeled by the 6-parameter analytical function [9], and the following six parameters were determined for each modeled twitch (Fig 2A): -the half-contraction time, the time from the start of the contraction until the moment when the twitch force reaches one half of its maximal value; T c (j) -the contraction time, the time from the start of the contraction until the moment when the twitch amplitude reaches its maximal value F max (j) (i); T hr (j) (i)-the half-relaxation time, the time between the start of the contraction and the moment when during the relaxation, the twitch force decreases to F max (j) (i)/2; T tw (j) (i)-the duration of the current contraction, i.e. the time from the moment when the contraction starts until the moment when the force decreases to 0.01% of F max (j) (i). The values of these parameters for the first decomposed twitch for the used in this study MUs are given in S1 File Table A. Besides these parameters, for each successive contraction, the level of the force at which it starts, F tetmin (j) (i) (Fig 2B), was also calculated from the force recording. The last parameter, F mftf (j) , was determined for each MU from the recorded tetanus evoked at 150 Hz as the maximum force of this tetanus and this value was accepted as the maximum tetanic force that a MU is capable of developing. F res is the force that the MU can still develop at its current state before reaching its maximum force of fused tetanus.
The lead time has a low variability [17], so for modeling purposes it was decided to keep this value unchanged, equal to the value obtained for the first decomposed contraction. Namely, it was accepted that T lead (1).  (1), for all 30 MUs. For each MU a linear approximation of the points (marked with different blue symbols for S MUs, different red symbols for FR MUs and different green symbols for FF MUs) was found, and the angles α (j) between the approximation lines and the Y axis were calculated. This angle was the smallest for the slow MUs, bigger for the FR MUs, and the largest (even over 90°) for the FF MUs. Values of this angle did not overlap for slow and fast MUs. From physiological point of view, the lines in the Fig 3 mean that for all slow and all FR MUs the amplitude of the next contraction within a tetanus increases when the developed force is closer to the maximal fused tetanic force, and this increase is higher for slow MUs with smaller angle α (j) . The amplitudes, however, have a limit, different for different MUs. This limit was not reached for slow MUs, probably because the tetanic curves were moderately fused and the capacity of slow MUs to develop more force is larger than that of fast MUs. Smaller force amplitude of a contraction for slow MUs is observed when the level of the force at which the contraction starts is far from the maximal possible force for a MU. For some FF MUs, for which α (j) >90 0 , the mentioned force dependencies are opposite-the amplitudes of the contractions increase with the increase of the residual force. The normalized residual force, however, cannot be very high. For the chosen set of FF Mus, the value of the coefficient F res (1) should not exceed 4. Fig 3 shows that the calculated angles α (j) are specific for each MU, meaning that the development of the tetanic contractions for each MU has its own course. Moreover, clear differences are visible between angles for slow and fast MUs. For these angles different dependencies were checked between different parameters of the decomposed twitches. The relationship between α (j) and the parameter F mftf (1) (Fig 4) was the most perspective for modeling. Moreover, the two parameters, F mftf (j) and F max (j) (1), i.e. the maximum possible MU force and the amplitude of the single twitch, were easily definable and specific for each MU. Several approximations of the above relationship with different linear and non-linear functions were tested using MATLAB functions. The power model shown in Fig 4 with square symbols was determined as the most suitable for our purposes since the calculated root mean square error was (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 F max (j) (i)/2; T tw (j) (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 F max (j) (i). The equation describing this bell-shape 6-parameters curve are given in [9]; B. Parameters measured for tetanic contractions presented on a part of the unfused tetanic curve (left) and the maximum fused tetanus (right). F mftf (j) -the maximal force that a MU develops during stimulation at 150 Hz stimulation frequency (in the fused tetanus). F tetmin (j) (i)-the force level at which the i-th contraction starts; F res doi:10.1371/journal.pone. 0162385.g002 (1), for all 30 MUs. 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-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-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 α (j) (j = 1, 2,. . ., 30) for each MU were calculated between these lines and the ordinate.
Having the specific angle α (j) for a MU, and constituting the equations of the straight lines shown in Fig 3, the amplitude of each i-th contraction within a tetanic curve could be calculated for each j-th MU by the following equation: The next two parameters of the successive twitch-like responses-contraction and halfrelaxation times-were hypothesized to depend on the actual MU force level at which the i-th contraction started, i.e. F tetmin (j) (i) (Fig 2B). The plots in Fig 5 suggest an approximation of the points with linear functions. Thus, the following equations for calculation of these two parameters could be written: The choice of the linear models shown in Fig 5 was confirmed by the high correlation coefficients, which were 0.8059 for T c and 0.7314 for T hr .
The remaining two parameters, i.e. the half-contraction time and the duration of the twitch, were supposed to be changed proportionally to the contraction and the half-relaxation times, respectively. Hence, the following equations were used for these two parameters: The tetanus force modeling procedure The Eqs (2), (3) and (4) include as a predictor the force level at which each contraction starts, i.e. F tetmin (j) (i). This is why first, the input parameters for modeling an unfused tetanic force curve for the j-th MU were specific for each MU, i.e. the force parameters F mftf (j) and F max (j) (1) and the twitch time parameters T lead (1), which could be estimated from the first decomposed contraction or from the recorded individual MU twitch, or were specific for the applied stimulation pattern, i.e. F tetmin (j) (i) (Fig 2). The latter parameter was at first calculated as local minima of the recorded tetanic curves. Having these parameters, the angle α (j) was calculated using the Eq (1). Then, the remaining parameters of the models of the successive contractions were calculated using Eqs (2)(3)(4)(5)(6). Note that the lead times were constants for all contractions of each MU. The twitches were modeled using the 6-parameter analytical function [9] and were summated according to the actual stimulation pattern. The obtained modeled force curve was compared with the respective experimental recording, and the two coefficients for estimation of the similarity between curves were calculated ( Table 1). The high accuracy of the modeling was demonstrated by the values of the two fit coefficients (Table 1, FitCo2 and AreaCo2).
However, the weakness of this modeling approach was the necessity of data concerning the actual force level, i.e. F tetmin (j) (i). For modeling purposes, this is a big obstacle when experimental data is missing for a given MU. That is why this parameter was also calculated in the next modeling step.
Instead of taking the parameters F tetmin (j) (i) from the recorded tetanic curves, F tetmin (j) (i) can also be predicted using the preceding modeled twitches (Fig 6). If the i-th twitch-like contraction of the j-th MU has the following force as function of time t (for t [0, t exp ], where t exp is the time duration of the experimental curve): where t IMP(n) is the time moment of the n-th pulse (n = 1, 2,. . ., 41) in the pattern used for the j-th MU. This mathematical process is visualized in Fig 6 where the first five contraction models are summed, the force developed by the MU at the moment of the 6-th pulse, F pr tetmin is calculated and used in Eqs (2), (3) and (4) as a predictor in order to predict the 6 th contraction.
For this type of more complex force prediction, the so called "full prediction", the procedure is the following: the α angle is calculated by the Eq (1), then the first contraction is modeled having F tetmin (j) (1) = 0; the force level F pr tetmin (j) (2) is computed from the first contraction as the force level at the moment when the second pulse is delivered, and this value is used instead of F tetmin (j) (2) in the Eqs (2)(3)(4). Hence, the second twitch model is generated. The first and the second models are added and the force level at which the third pulse comes is calculated. This is the value F pr tetmin (j) (3), which is used instead of F tetmin (j) (3) in the Eqs (2)(3)(4). This procedure continues until the last, 41 st pulse. This way a tetanic force can be predicted from 6 available parameters of the single twitch, the maximal fused tetanus force and a given stimulation pattern.

Results
The curves obtained by summation of the models of all 41 decomposed contractions well resembled the recorded tetanic curves for all 30 MUs. The mean value of the fit coefficients between the recorded and the modeled curves was 98.6829 (range 98.0604-99.2647). The mean value of the AreaCo was 1.0068 (range 1.0003-1.0154), i.e. the modeling error was below 2%. There were no differences in the preciseness of the modeling between different types of MUs. This confirmed the preciseness of the decomposition of all 30 tetanic curves and completeness of the teaching database.
When using all input parameters (S1 File Table A) for the generation of the prediction model curves, including the force level at which each successive contraction starts taken from (1) and F max (j) (1), respectively. The equation of the straight line is: y = p 1 +p 2 x, and p 1 = 2.397, p 2 = 0.3509; blue circles-S MUs; red circles-FR MUs; green circles-FF MUs. The two twitch time parameters used to calculate the data presented on the ordinate are illustrated in frames left to the axis on the example of a series of decomposed twitches (red lines indicate the contraction time for the first and the i-th twitch in A, and the half-relaxation time for the first and for the i-th twitch in B). 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 recording (red lines indicate amplitudes of the first twitch and the force level at which the i-th contraction starts).
doi:10.1371/journal.pone.0162385.g005 the force recordings, the tetanic curves modeled by using the approximation approach, i.e. applying the Eqs 1-6, were much closer to the recorded ones than to the curves obtained by summing equal twitches for all MUs (Fig 7). Note that the IPIs for the stimulation patterns used for these three tetanic curves are given in S1 File Table B. The visual observations from these three plots were confirmed by the values of the two coefficients given in Table 1. The number 1 was assigned to the coefficients estimating the similarity between the recorded curve and the curve obtained by summation of equal twitches. The number 2 was assigned to the coefficients, estimating the similarity between the recorded curve and the curve obtained by Table 1. Calculated angles and coefficients reflecting the similarity between the experimental curves and the predicted curves for all 33 MUs. FitCo1 and AreaCo1 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. α1 (j) is the angle calculated by using the Eq (1). FitCo2 and AreaCo2 are the coefficients for the experimental and modeled curve (the calculated experimental values of F tetmin (j) (i) are used as input parameters for the prediction). FitCo3 and AreaCo3 are the coefficients for the experimental and modeled curve using the same angles α1 (j) but for full prediction algorithm (i.e. values of F tetmin (j) (i) are also predicted consecutively). α2 (j) is the angle obtained by a sensitivity analysis so that the modeled curve is the most similar to the experimental one, and FitCo4 and AreaCo4 are the respective coefficients. α3 (j) is the improved angle using the Eq (11) (see Fig 8), and FitCo5 and AreaCo5 are the respective coefficients.

MU FitCo1
AreaCo1 α1  the summation of the contraction models, whose parameters were predicted by the approximation. FitCo1 was always lower than FitCo2 while AreaCo1 was always bigger than AreaCo2 (with the exception of FF8), and AreaCo1 was usually closer to 1. The most apparent differences were observed for slow MUs. For slow MUs the difference between the two fit  Comparison between the recorded tetanic curves (red), the force curves obtained by summation of equal twitches according to the same stimulation pattern (green), and the curves predicted by the approximation approach (blue) using the angles α1 (i) from coefficients, FitCo1-FitCo2, varied between 21.4256 (for S7) and 51.2561 (for S4, this was also the maximum value for all MUs). Smaller values for FitCo1-FitCo2 were observed for FF MUs; the minimum value was 0.4302 (for FF8), and the maximum value was 23.1124 (for FF3). The minimum absolute value of AreaCo1-AreaCo2 was 0.0522 (for FF8), while the maximum was 6.3325 (for S6). To summarize, the force curve predicted by the model using the new proposed approach was always considerably closer to the recorded tetanic force than the one obtained by summation of equal twitches. This conclusion was more evident for the slow MUs. Table 1 and Fig 3 show the range of the angle α (j) , which varied from 53.8892°to 95.8677°, and increased from slow MUs to fast ones. For three FF MUs, this angle was even bigger than 90°. When this angle was equal to 90°the approximation line for a particular MU was horizontal and all successive contractions had one and the same force amplitude, independently of the stimulation pattern. These twitch-like contractions only differed in their shapes by a prolonged duration, because of the dependencies shown in Fig 5 and the Eqs (3-6). Within the population of the 30 MUs taken as input data, the relationship between the maximal possible tetanus force of a MU and the amplitude of its single twitch force (hence, F mftf (1)-see S1 File Table A), was from 1.6260 (for FF6) to 14.8666 (for S1), and this relationship explicitly determines the angle α (j) (see Eq (1)) for a particular MU.
In the majority of FF MUs, some values of F max (1) were below 1 (Fig 3). This means that some decomposed contractions had amplitudes lower than the first contraction. This never occurred for slow and FR MUs.
The plots in Fig 3 depend on two important parameters, specific for each MU-the maximal possible tetanus force that a MU can develop and the amplitude of the first decomposed twitch. Both are subjects of error during modeling. This especially refers to changes of the parameters of the twitches due to physiological processes, such as a potentiation at the beginning of an experiment or fatigue at the end. As an effect, the recorded force can vary within the experiment. On the other hand, the maximal tetanic force can decrease progressively during an experiment due to the development of fatigue, especially in FF MUs.
To check whether the proposed mathematical approach could be applied for various MUs in rat medial gastrocnemius (i.e. MUs not included in the initial database) three additional MUs, one of each type (S11, FR11 and FF11 in Table 1 and in S1 File Table A), previously not taken for decomposition and modeling procedures, were chosen for verification. The recorded tetanic curves of these MUs had similar fusion indices, but stimulation patterns (given in S1 File Table B) were different from all those previously used. The predicted force curves closely resembled the recorded curves for S11, FR11, and FF11 MUs. This was confirmed by the coefficients FitCo2 and AreaCo2 given in Table 1.
The next step in the modeling was an attempt to fully predict the course of a MU tetanic contraction, calculating the consecutive force levels at which the subsequent pulses were delivered, using Eq (8), but not taking them from the recorded tetani. Hence, only the parameters specific for a given MU and the stimulation pattern would be necessary for a tetanic force prediction. This so called "full prediction" was made for all 33 MUs. FitCo3 and AreaCo3 in Table 1 reflect the similarity between the recorded force curve and the one modeled by the full prediction algorithm. The recorded and the predicted curves for the S11, FR11 and FF11 MUs, which are outside the database, are shown in Fig 8. Comparing FitCo3 and AreaCo3 with FitCo2 and AreaCo2 from Table 1 it can be concluded that the preciseness of the prediction changes for the worse since additional inaccuracies were put into the calculations-the force level at which the next contractions start. Possible mistakes could be also due to the calculation of the angle α (j) from Eq (1), since this parameter considerably influences the force amplitude of the consecutive contractions. Small changes in this angle could potentialy increase or decrease the tetanic force. Moreover, the dependence shown in  Table 1, the angle α1 (j) , calculated using Eq (1), was changed within a suitable limit, and a value giving the best prediction was chosen as the angle α2 (j) . The limits for the sensitivity analysis were the integer value of α1 (j) plus/minus 20°with steps of 0.1°. The best prediction was chosen on a basis of both coefficients FitCo4 and AreaCo4 -they had to be as close as possible to 100% and 1, respectively. Table 1 shows values of α angle corresponding to the most accurate prediction of the tetanic curves. For all slow MUs (except S11) this angle was decreased, and for one MU it was decreased even by 12°. For nearly all fast MU (except for FR9, FF2, FF6, FF9 and FF11) this angle was increased, and for FF4 it was increased even by 19°.
These observations were used to improve the Eq (1) . Fig 9, where these angles α (j) , obtained after the sensitivity analysis, are plotted versus F mftf , shows a better fit function. An ameliorated equation for the dependence in Fig 9 was The full prediction was performed again for all 33 MUs, using the new equation. The final values of the angle were calculated using the Eq (9), i.e. α3 (j) in Table 1. The coefficients FitCo5 and AreaCo5 are in general better than FitCo3 and AreaCo3, although there are not as good as FitCo2 and AreaCo2. This allows us to suppose that the precise values of F tetmin (j) (i) are very important for predicting the successive contractions and that a further sensitivity analysis with respect to the contraction and the half-relaxation times is necessary.

Discussion
This paper has two achievements. First, the proposed approach can be used to develop more realistic muscle models composed of a set of MUs, which may be used to study processes of control of the skeletal muscle force. Second, the differences in the force development during repetitive activation of MUs from various types have become clearer, and the separation between the three types of MUs can now be made based on the new observations, namely the  values of the angle α (j) (Figs 3 and 9). The application of the angle α is the major novelty of this muscle model. It appears to be a parameter specific for respective types of MUs and enables us to predict the course of unfused tetanic contractions using the same equations and MU quantities for arbitrary stimulation patterns with variable interpulse intervals. Fig 3 shows that each MU has its own specific angle. The values for α1 (j) for slow MUs ranged from 53.8892°to 62.2282°(S11 seems an exception), and for fast units-from 66.2353°to 95.8677°, so these ranges were not overlapping ( Table 1). The values of the angle α1 (j) overlapped for FF and FR MUs, but on average, were higher for FF units. The same conclusions can be made for α2 (j) and α3 (j) . For slow MUs, α2 (j) is between 44.5°and 60.1°, and for fast MUs, α2 (j) is between 74°and 96°. For slow MUs, α3 (j) is between 50.1630°and 60.2272°, and for fast MUs, α3 (j) is between 74.3563°and 100.5887°. In general, this angle determines whether the amplitudes of the successive contractions within an unfused tetanus would increase (α (j) <90°), remain constant (α (j) = 90°) or decrease (α (j) >90°) during tetanic force development. The range of the values of this angle depends on the relationship F mftf (j) /F max (j) (1), so one may expect that for another population of MUs (other muscles, muscles in various species), ranges that are different from those presented in Table 1 would be obtained.
The MU force varies within a range from the single twitch amplitude up to the maximum (fused) tetanus force. This range is usually expressed by a twitch-to-tetanus force ratio, and this ratio is lowest for S MUs and highest for FF units. For the rat medial gastrocnemius, the mean values of this ratio amount to 0.13±0.05, 0.19±0.06 and 0.28±0.08 for S, FR and FF MUs, respectively [18]. Accepting that the minimal value of the twitch-to-tetanus ratio is about 0.05, and that p denotes this ratio, the inequality 0.05 p 1 can be written since the twitch amplitude has to be lower than the maximal possible MU force. This can be rewritten as 0.05 F max (j) (1)/F mftf (j) 1. Then because of the improved Eq (9), a ðjÞ ¼ 117:2ð1=pÞ ðÀ 0:3144Þ : Hence, the maximal value of the angle is 117.2 0 and the minimal one is 45.6966 0 . So, the Eq (9) has clear physiological basis.
Neither Eq (1) nor Eq (9) imposed a mathematical limit for the angle. However, since an MU force may vary from a minimum activity during a single twitch up to the maximum force in the fused tetanic contraction, a limit based on the maximum possible tetanic force should be taken into consideration during modeling more fused tetani. From a physiological point of view, the power model lines in Fig 4 and Fig 9 have to be limited by vertical lines crossing 1 at the abscissa, since the amplitudes of the individual twitches cannot exceed the maximal tetanus force.
Nevertheless, the proposed equations predict some limits in the force development. They have physiological interpretation, too. For example the Eq (2)  (i)>0, i.e. in case when the next pulse comes before the full relaxation of a MU, then the amplitudes of the next contractions would always be bigger than F max (j) (1), since cot(α (j) )>0. If α (j) = 90°, the amplitudes of the next contractions would always be the same as F max (j) (1). The third case, when α (j) >90°, seems especially interesting, since the amplitudes would decrease. Three examples are given in Fig 10 for the three border variants of the angle: α (j) = 45 0 , α (j) = 90 0 and α (j) = 117.2°. As can be concluded from Fig 3 and Table 1, for all slow and FR MUs, α (j) <90 0 . In such a case ( Fig 10A) the amplitude of a contraction within a tetanus increases when the level at which the contraction starts is higher. If F tetmin (j) (i) = 0, i.e. in the case when a MU is fully relaxed, the contraction will have the same amplitude as the single twitch (the first point in Fig 10A). Upper limit of F max (j) (i) and F tetmin (j) (i) has to be imposed further by inequality constraints for more fused tetani. When α (j) = 90 0 ( Fig 10B) the amplitudes of all contractions of this MU will be equal, i.e. F max  (1) when the next stimulus is delivered, the MU will be not able to respond with a new contraction, i.e. the MU will not generate more force.
The MU force can be regulated by changes in the motoneuronal firing rate, and during voluntary activity, motoneurons generate trains of pulses at variable intervals, which substantially influence the development and profiles of tetanic forces [19][20][21]. In our study, 13 different patterns with different mean frequencies were used for stimulation of the MUs. This diversity of the stimulation patterns was applied in the experiments to ensure that statistical data is reliable for approximation (S1 File Table B). Moreover, the approach was validated using the tetanic force curves obtained by stimulations with three additional patterns, different from those used previously. It has to be mentioned also that the range of the contractile parameters of the 33 investigated MUs (S1 File Table A) covers the range reported for MUs in rat medial gastrocnemius muscle with respect to all basic parameters [22].
The derived approximation Eqs (1-4) are based on 30 decomposed tetanic curves (Figs 3, 4 and 5 for dependencies). The most critical appears to be the relationship shown in Fig 4. The possible errors come from two sources: (1) the amplitudes of the single twitch and that of the first decomposed contraction are different for a MU depending on fatigue, potentiation, previous stimulations, etc. [23][24][25]; (2) the maximal tetanic force, especially for FF type MUs, may decrease due to fatigue and the estimation of this maximal force may also evoke an error. This explains why other angles-α2 (j) -can significantly improve the prediction, and this was the reason for proposing the new, improved, power model by Eq (9).
Another reason for the impreciseness of the full prediction, using the calculated angle α1 (j) ( Table 1), is that the predicted values of F pr tetmin (j) (i) are slightly different from the actual force levels at which the successive contractions start. Each of them starts with a delay of T lead (j) (i) (Fig 6). This lead time is not known in advance. It can be considered for the further improvement of the algorithm taking its value for the first contraction, i.e. T lead (j) (1), which can be added to t IPI(n) in Eq (8).
The input data used for approximations were the results of decompositions of tetanic contractions recorded for a group of three types of MUs, for different stimulation patterns, with different mean frequencies, but at similar fusion indices and similar relative force levels. An open question remains: whether the approach would be suitable for force curves with considerably different fusion indices (for example, for nearly fused tetanic curves) or for stimulation with equal IPIs. As discussed above, the derived equations might be not applicable for MUs of other rat muscles or for muscles in other species. Nevertheless, one may expect that the equations for force modeling in other muscles would be based on a similar main idea when matched to specific databases for these muscles. Irrespective of the above restrictions, the proposed new approach for predicting the successive twitch-like responses fusing into tetanic contraction has a general character, due to the fact that the same algorithm can be used for all MU types within a muscle. It is also important that, with a progress in motor units' research, the described model could be further suitably expanded by adding new modules related to the non-linearity in the summation of forces of many MUs, to potentiation observed in fast MUs during their long-lasting activity and to fatigue in fast MUs. At present, there is lack of sufficient data concerning these phenomena.
Supporting Information S1 File. The six parameters of all motor units used in this study allow to graphically present each individual (or first decomposed) twitch using the analytical function given in Raikova et al. (2008). They constitute the database for the prediction algorithm presented in the paper (S1 File Table A). S1 File Table B presents examples of the used stimulation patterns, i.e. the IPIs between the subsequent 41 pulses. S1 File Table A. The basic parameters used when the approximation was built and the reconstruction was made (for explanation, see Fig 2). Consecutive rows present data for the 10 slow MUs, the 10 FR MUs and the 10 FF MUs used as an input database. The last three rows show the data for the three additional MUs of three types (S, FR and FF), which were used for verification of the approach. The 10 MUs of each type are listed in order of increasing forces of the first twitch, F max (j) (1). S1 File Table B. Interpulse intervals for stimulation patterns used for experimental tetani shown in Fig 7 and