Effect of synchronization of firings of different motor unit types on the force variability in a model of the rat medial gastrocnemius muscle

The synchronized firings of active motor units (MUs) increase the oscillations of muscle force, observed as physiological tremor. This study aimed to investigate the effects of synchronizing the firings within three types of MUs (slow—S, fast resistant to fatigue–FR, and fast fatigable–FF) on the muscle force production using a mathematical model of the rat medial gastrocnemius muscle. The model was designed based on the actual proportion and physiological properties of MUs and motoneurons innervating the muscle. The isometric muscle and MU forces were simulated by a model predicting non-synchronized firing of a pool of 57 MUs (including 8 S, 23 FR, and 26 FF) to ascertain a maximum excitatory signal when all MUs were recruited into the contraction. The mean firing frequency of each MU depended upon the twitch contraction time, whereas the recruitment order was determined according to increasing forces (the size principle). The synchronization of firings of individual MUs was simulated using four different modes and inducing the synchronization of firings within three time windows (± 2, ± 4, and ± 6 ms) for four different combinations of MUs. The synchronization was estimated using two parameters, the correlation coefficient and the cross-interval synchronization index. The four scenarios of synchronization increased the values of the root-mean-square, range, and maximum force in correlation with the increase of the time window. Greater synchronization index values resulted in higher root-mean-square, range, and maximum of force outcomes for all MU types as well as for the whole muscle output; however, the mean spectral frequency of the forces decreased, whereas the mean force remained nearly unchanged. The range of variability and the root-mean-square of forces were higher for fast MUs than for slow MUs; meanwhile, the relative values of these parameters were highest for slow MUs, indicating their important contribution to muscle tremor, especially during weak contractions.


Introduction
Most studies of motor unit (MU) firings have revealed the existence of a certain level of synchronization between the firings of motoneurons innervating the same muscle [1][2][3][4]. Two concepts for long-and short-term synchronization can be found in the literature. Long-term synchronization with greater latencies beyond ± 20 ms was reported by Datta [1,[4][5][6][7]. The possible mechanism of this kind of synchronization could be explained as interactions occurring between the stretch reflex loop and the recurrent inhibition. Long-term synchronization has been reported to be relatively rare in comparison to short-term synchronization [4], which was reported to be a peak in the cross-interval histogram centered about a zero-time delay (0.5 ± 2.9 ms). Short-term synchronization is attributed to last-order projections that provide common, nearly simultaneous, excitatory synaptic input across motoneurons [3,8,9,10], generating a narrow peak around the origin of the cross-correlogram of MU discharges [1,8,11,12]. Therefore, the narrow synchronous peak principally reflects shared, monosynaptic projections to motor neurons from corticomotoneuronal cells via the lateral corticospinal tract [13].
In humans, the MU synchronization was shown to be stronger during voluntary muscle activation than during reflex activation [14]. At the same time, synchronization tends to be higher in more distally located muscles, while the greatest synchrony has been most often found in the intrinsic muscles of the foot rather than in the hand muscles [3,15]. However, the level of synchronization between MUs could be influenced by numerous factors, such as the examined task, the muscles involved in the task, and the type of habitual physical activity performed by the individual [6][7][16][17][18][19]. For example, the level of synchronization was reduced between MUs in the hand muscles of individuals who required greater independent control of the fingers. This included musicians [18] and the dominant hands of control subjects [7]. On the other hand, MU synchronization was found to be greater in the hand muscles of individuals who consistently performed strength training [18,20] or during tasks that demanded attention [21]. The enhancement of MU synchronization was observed after daily exercise involving brief periods of maximal muscular contraction [20] and contributed to traininginduced increments in muscle strength [22]. Greater synchronization has also been noted in fatigued muscles [23]. Reports regarding the relationship between physiological tremor and synchronization are inconsistent: most of them have linked tremor with an increased level of synchronization [23][24][25][26], while others have suggested no significant associations between the tremor amplitude and the level of MU synchronization exist [18].
It has been assumed that muscle can produce smooth contractions due to asynchronous discharges of motor neurons [24]. Yao et al. [22] revealed that MU synchronization increased the variability in the simulated force but not the average force. Synchronization was also shown to increase the estimated twitch force and decrease the contraction time of the MUs [27], though it was likely a result of inaccurate measure of MU twitch properties in that model.
In the majority of skeletal muscles, three types of MUs have been distinguished and their contractile properties, including the force-frequency of stimulation relationship [28] and sensibility to changes in stimulation pattern [29,30], were found to vary considerably. In several studies, the effects of the synchronization of MU firings were modeled [22,31,32]; however, these models did not analyze the specific effects attributable to different types of MUs. In our previous paper [33], a model of the rat medial gastrocnemius muscle consisting of 30 MUs [10 MUs each of the fast fatigable (FF), fast resistant to fatigue (FR), and slow (S) types] was proposed and the effects of synchronous and asynchronous stimulation of MUs were investigated. It was concluded that the activation of MUs at variable interpulse intervals, delivered to each MU asynchronously, resulted in smaller force oscillations. However, the study did not assess the effects of synchronization between pairs of individual MUs nor the effects of the synchronization of three types of MUs.
A recent model of the rat medial gastrocnemius muscle [34] provided methodology by which to identify the role of each of three MU types (FF, FR, and S) in the production of muscle force. In the present study, the same model was adopted as a tool for simulation of four modes and three time levels of synchronization. The aim of this research was to reveal the effects of synchronization of firing within the three types of MUs on the force variability and the force mean spectral frequency and to compare these effects between different types of MUs and the whole muscle. The implication of the results for explanation of tremor at various levels of the muscle force was discussed.

Muscle model
This study applied a model of the rat muscle gastrocnemius based on excitability and firing frequencies of motoneurons, contractile properties, and the number and proportion of MUs in the muscle [34]. Briefly, the model consisted of 57 MUs, including eight S, 23 FR, and 26 FF MUs, respectively. Their basic contractile properties are presented in Table 1. As input data, this set of MUs, recorded in physiological experiments, was selected and their twitches were precise modeled by a six-parameter analytical function dependent on time [35]. These twitches are given in Fig 1 in Raikova et al. (2018) [34]. The muscle force was calculated as the sum of forces of all active MUs and the process of force regulation was set according to the commondrive hypothesis [36]. The muscle unfused tetanus was calculated following the application of a train of irregular stimuli and was simulated using an analytical approach described in previous researches [34,37]. The tetanic force of each MU was calculated as a sum of unequal twitches analytically described by the 6-parameters analytical function. These parameters were changeable and were computed using regression equations derived by decomposition of tetanic curves recorded in physiological experiments for 30 MUs [37]. Meanwhile, the scheme of MU firing was adopted from Fuglevand et al. [31]. Table 1. Main input data for the muscle model: values of the basic contractile properties of the 57 MUs (slow, S1-S8; fast resistant to fatigue FR1-FR23; fast fatigable FF1-FF26) selected for the muscle modeling and their firing properties. T c -the contraction time; T hr -the half-relaxation time; T tw -the duration of the twitch; F max -the maximum force of the twitch; F mftf -the maximum force of the fused tetanus; meanfr-the mean frequency of rhythmic firing; minfr-the minimum rhythmic firing frequency; maxfr-the maximum rhythmic firing frequency. Data are taken from Table 1 in Raikova et al. [33]. In the present study, the excitation signal was simulated ( Fig 1A) as consisting of two smooth logarithmic parts during the rising and falling slopes of the muscle force (each lasting 1000 ms) and a straight line during the steady state of the muscle (lasting 2000 ms). The shape of the signal waveform was designed to better approximate more realistically a course of excitation input to motoneurons, avoiding sudden changes occurring in any trapezoidal signal used previously. This study considered only one excitation level, corresponding to 100% of the activation signal, ensuring that all MUs were activated during the steady state of the muscle to enable a thorough analysis of their synchronization. The program for simulation of the force MUs and the muscle force accepted the same MU firing frequencies as previously described [34] ( Table 1). The interpulse intervals (IPIs) were calculated by generator of uniformly distributed pseudorandom numbers within the interval of ±4 ms around the IPI corresponding to the mean frequency of each MU. The simulated train of firings at random IPIs are presented in a supplementary table (S1 Data). IPIs histograms are also presented in a supplementary figure (S1 Fig) demonstrating IPIs uniform distribution. Finally, the model generated the output forces for different MUs (S, FR, and FF) and the whole muscle, as illustrated in Fig 1B (sampling frequency fs = 1 kHz) and raw data presented in the supplement (S2 Data). This was further denoted as the basic (non-synchronized; NS) model, to which no attempts of manual changes of MU firing for synchronization were applied. The force signals were analyzed during the steady-state periods (2000-4000 ms). Their power frequency spectra were calculated by using fast Fourier transform (FFT) over nf = 2048 points, thus achieving a spectral resolution Δf = fs/nf = 0.49 Hz (Fig 1C).

Simulation of MU synchronization
The NS firings of all 57 MUs during the muscle steady state are shown in Fig 2 and presented in a supplementary table (S1 Data). These patterns were further modified to simulate different types and levels of synchronization. The synchronization was applied to a specific pair of MUs (named MU1 and MU2) so that the pulses of MU2, which fall within a predefined time window, Δt, around the pulses of MU1, were changed to coincide with those of MU1. Three time windows with Δt = ± 2, ± 4 or ± 6 ms were used to simulate three levels of MU synchronization, mimicking weak, modest, and strong synchronization, respectively. The synchronization scheme is illustrated in Fig 3, showing that the larger the time window was, the greater number of MU pulses were shifted to and synchronized with the reference MU.
Four methods of synchronization (Methods 1-4) were applied. In all methods, the synchronized MU pairs were chosen only encompassing the same physiological type (S and S, FR and FR, FF and FF), i.e., synchronization was not induced between MUs of different types.

Estimation of MU synchronization
Temporal correlation of MU pulses. The MU pulses were represented as an MU binary (MUB) sample series with a constant sampling period of 1 ms and binary amplitude of 0 or 1, where "0" indicated a non-active state and "1" indicated the presence of a pulse-active state. The duration of the pulse-active state was set to 1 ms, overlaying one sampling period. MUB series were represented with a total of 2000 samples during the steady state of the muscle from 2000 ms to 4000 ms, as depicted in Figs 1 and 2 for the MUs in the basic, NS model.
The temporal correlation between the binary sample series of two MUs (MUB1 and MUB2) was computed with the normalized Pearson's correlation coefficient ranged in the interval 0% to 100%, according to the following formula: ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi X where i denotes the sample index of the MUB series, considering a sampling period of 1 ms. The correlation coefficient (corMU) is a standard measure of similarity between sample series data in the time domain. Transferring this knowledge to the MUB data, corMU is representative of the temporal synchronization of two MU firings such that 100% corresponds to a complete coincidence between all firing pulses in MU1 and MU2 and 0% corresponds to no coincidence between any firing pulse in MU1 and MU2. The normalized value of corMU does not depend upon the length of the estimated MUB time series, the number of firing pulses, or the mean firing rate. This is an important benefit of the normalization, which would prevent bias in corMU estimation, considering that MUs in different physiological types have different mean firing rates.
Cross-interval synchronization index. The synchronization between the firing patterns of two MUs (MU1 and MU2) was estimated by an analysis of their cross-intervals using CI x (MU1,MU2) = {t1 x − t2 xy } computed as a pair-wise difference between the times of occurrence of all reference MU1 pulses t1 x = {t1 1 ,t1 2 ,. . ..t1 nMU1 } and their corresponding closest neighbors among MU2 pulses t2 xy 2t2 y = {t2 1 ,t2 2 ,. . ..t2 nMU2 }. The latter were found by the minimization criterion t2 xy ¼ arg min y¼1;2;::nMU2 fjt1 x À t2 y jg and respected the condition to overlay only firings during the steady state of the muscle, i.e., t1 x ,t2 y 2[2000ms; 4000ms]. By definition, the CI x (MU1,MU2) vector length was equal to the number of pulses in the reference MU (nMU1). CI values could be negative, zero, or positive when an MU1 pulse was respectively preceding, coinciding with or following its neighbor MU2 pulse, as illustrated in Fig 4. The distribution of cross-interval values of two MUs, CI(MU1,MU2), was estimated by means of a cross-interval histogram with a bin-width resolution of 1 ms and bin centers in the range of ± 15 ms. The bin values represented the relative probability (p bi ) of having a CI observation within a specific bin interval (bi): accepting the sum of all bin values equal to 1: where c bi is the count of CI (MU1,MU2) values in bin bi and the denominator is the number of elements in the input data, equal to the number of reference MU pulses (nMU1). Derived from the cross-interval histogram, the synchronization between the firing patterns of MU1 and MU2 was estimated by the relative probability p b0 in the central bin (b0 = ± 0.5 ms), equivalent to the relative frequency of coincidence between MU1 and MU2 pulses related to the reference number of pulses: Given a total number of N = 57 MUs, there could be derived a total of N-1 cross-interval vectors CI(MUi, MUj) for any given pair of MUs, where i, j = 1, 2,.., N, and i6 ¼j. Further, a cross-interval synchronization index (CISI) was defined for each reference MUi pattern to accumulate the relative probability of pulse coincidences in all MUi pairs (N−1) observed in the central bin: CISI has a normalized value from 0 to 100% with 0% corresponding to no coincidence and 100% corresponding to a complete coincidence between the patterns of the reference MUi and all other paired MUs. The adopted CISI normalization to both number of MU pairs (N) and number of reference firings (nMUi) was implemented to reject the influence of the size and type of the studied MU population. Force parameters. Standard force parameters of the different MU groups (S, FR, and FF) and the cumulative force of the whole muscle ( Fig 1B) were calculated as follows: • Force root-mean-square (RMS) level: rmsF ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi X n i¼1 where F i denoted the force signal samples, taken with a sampling period of 1 ms during the steady state of the muscle from 2000 ms to 4000 ms, including a total number of n = 2000 samples.
Furthermore, a variance based parameter, namely percentage of Variance Accounted For (VAF), was computed to account the effect of MU synchronization on the variance of differences between synchronized forces (FS) and the non-synchronized forces (FNS) of the basic, NS model: where i denotes the force signal samples from 2000 ms to 4000 ms, var(FS i -FNS i ) denotes the total sum of squared differences between forces and var(FNS i ) denotes the total sum of squares taken with respect to the mean value. VAF was used in literature to estimate the synchronization between forces of different muscles [38,39,40]. VAF of two equal forces is 100%. VAF is expected to decrease if the variance of the differences between two forces grows due to synchronization. Additionally, the force power spectral density (PSD) of different MU groups (S, FR, and FF) and the cumulative muscle force (Fig 1C) was used for the calculation of the mean spectral frequency as follows: where nf is number of frequency bins in the spectrum (nf = 2048 as defined earlier), fi is the frequency of the spectrum at bin i of nf, and PSDi is the amplitude of the PSD at bin i of nf.

Weak synchronization of MU firings in the basic muscle model
The level of synchronization between MU firing patterns of the simulated basic rat muscle gastrocnemius model with 100% excitation containment and 57 MUs, over a two-second time period during the muscle steady state, was estimated by the two different synchronization indices in Table 2 (top row) and discussed as follows. First, corMU = 6.1% ± 2.8% (mean value ± standard deviation) shows a weak temporal correlation between the firing pulses of all MUs within the muscle, which was found to be lowest for MUs of type S (4.5% ± 2%) and highest for those of type FR (7.4% ± 2.9%). A comprehensive proof for the absence of noteworthy clusters with significant correlation between MUs of a specific type is illustrated in the corMU colormap in Fig 5A. Here, a random distribution of corMU values can be noted, overlaying the dark-blue colored area of very low pairwise correlations between 57 × 57 MUs, distributed on the x-and y-axes. The entries in the main diagonal should be ignored because each represents a MU compared with itself (corMU = 100%).
Second, CISI = 6.2% ± 0.4% (mean value ± standard deviation) suggests weak cross-interval synchronization between the firing pulses of all MUs within the muscle, without essential differences between MUs of different physiological types [the CISI mean value varied from 5.8% (S MUs) to 6.2% (FR and FF MUs)]. Evidence for missing synchronization between MU firings can be observed in the cross-interval histograms in Fig 6A, having a flat (uniform) distribution in the range of bin intervals [−6 ms; +6 ms] for all 57 MUs. Therefore, cross-intervals between firing patterns were equally probable within this bin range and no evidence for synchronous peaks could be identified in the case of any MU.

Stronger synchronization of MU firings in different synchronization scenarios
The aforementioned 57 MU firing patterns of the basic muscle model were modified according to 12 synchronization scenarios, i.e., four synchronization concepts (Methods 1-4) each applied within three synchronization time intervals (Δt = ± 2, ± 4, and ± 6 ms). The resultant average levels of synchronization between patterns of MUs of the same physiological type and within the whole muscle are estimated in Table 2. In all cases, certain increments of both indices for the level of MU synchronization (corMU and CISI) were assessed in comparison with their estimation for the basic NS model in the first row of Table 2. Therefore, it may be concluded that the simulation design achieved the general goal for inducing stronger synchronization between MU firings. More details on the observed MU synchronizations related to the computation of corMU and CISI are presented below.
• corMU: Different effects of the synchronization induced by Methods 1 through 4 could be tracked well on the corMU color map (Fig 5B, 5C, 5D and 5E), seen as clusters with strong correlations (corMU is from 30% to 100%). These clusters have different two-dimensional space distributions of the entries with maximal correlation, corresponding to the different concepts for synchronization of MU pairs in Methods 1 through 4, as follows: � Method 1: The synchronization between neighboring MUs is seen in Fig 5B as (Table 2). This result yields an increment from 37% to 67% of the correlations within MU groups relative to the basic NS model and can be denoted as the maximal synchronization level simulated in this study.
• CISI: The effect of synchronization induced by Methods 1 through 4 could be identified in the cross-interval histograms (Fig 6B, 6C, 6D and 6E) by the prominent peak in the central bin (± 0.5 ms). The larger is amplitude deviation from the uniform distribution in other bins, the higher the probability for synchronization of the respective MU to the firing pulses of other MUs. Different synchronization methods produce different amplitudes in the central bin, estimated by CISI in Table 2. Comparing the CISI values of all methods estimated with maximal synchronization interval Δt = ± 6 ms, we could deduce the following: � The lowest CISI mean value was found for S MUs (from 8.2% in Method 3 to 10.3% in Method 2, with the latter being up to 4.5% above the basic NS model).  Fig 7A was designed to show the effect of widening the time window for synchronization (Δt = ± 2, ± 4, and ± 6 ms) on the relative CISI change (ratio of synchronized vs. NS value). It shows that, generally, Δt = ± 2 ms leads to weak synchronization and slight increases in CISI by about 1.1 times (Methods 1-3) and 1.8 times (Method 4); Δt = ± 4 ms lead to maximal synchronization that is still less than two times (Methods 1-3) but about three times (Method 4); and Δt = ± 6 ms produced the maximal synchronization with notable CISI increment increases by up to three times (Methods 1 and 2) and up to five times (Method 4).

Maximal effect of MU synchronization on the force parameters
The forces produced by the muscle and different MU types before and after the application of different synchronization scenarios were estimated for a two-second period during the muscle steady state and the defined six basic force parameters (meanF, rmsF, rangeF, maxF, VAF and meanfreq) are presented in Table 3. For comprehension purposes, the representation of those parameters on the force signals and their PSD (power spectral density), is additionally illustrated in Fig 8. The comparison of the NS excitation to those achieved with different synchronization methods (Methods 1-4) is presented below.
• Force mean value: The synchronization had no effect on meanF parameter, which accounted the value of the mean force. It showed a negligible change (� 12 mN) before and, after the synchronization was applied, i.e., for the muscle force, meanF was varying from 4052 mN (NS) to a maximum of 4064 mN (Method 4, Δt = ± 6 ms) ( Table 3). This can be also tracked in Fig 8A, which presents no visible difference in the baseline value (red solid horizontal line) when comparing all forces placed in a row.
• Force RMS value: The synchronization had an important effect on the force variance, increasing the rmsF value by more than 50 mN. It could become as high as 129 mN for the muscle force (Method 4, Δt = ± 6 ms), considering its baseline NS value of 73.5 mN (Table 3). Additionally, Fig 7B is provided to show the relative rmsF change as a ratio of synchronized vs. NS value. Specifically, it shows that the maximal rmsF increment (about two times) could be achieved for the forces of two types of MUs (S, FR) following synchronization with Methods 1, 2, and 4, Δt = ± 6 ms. Considering the whole muscle, the observed maximal increment of rmsF was about 1.8 times, achieved using Method 4, Δt = ± 6 ms.
• Force min-max range: The synchronization had an important effect on the force variance, increasing the rangeF value by about 450 mN. It grew from 405 mN (NS) up to 850 mN for the muscle force after synchronization with Method 4, Δt = ± 6 ms ( Table 3). The rangeF ratio (synchronized vs. NS value) in Fig 7C shows that the largest rangeF increment (two to 2.6 times) was achieved for the forces of two types of MUs (S, FR) after synchronization with Methods 1, 2 and 4, Δt = ± 6 ms. Considering the whole muscle, the observed maximal increment of rangeF (about 2.1 times) was with Method 4, Δt = ± 6 ms. Although the observations concerning rangeF are similar to those of rmsF as was noted above, the relative and absolute  changes in rangeF values as an effect of synchronization were larger. This could also be visually confirmed by the force signals in Fig 8A (blue dotted lines show larger span than red dotted lines after synchronization, comparing all forces placed in a row).
• Maximal force: The synchronization had an important effect on increasing the maxF value by more than 205 mN, which could raise it from 4234 mN (NS) up to 4440 mN for the muscle force after synchronization with Method 4, Δt = ± 6 ms ( Table 3). The maxF ratio (synchronized vs. NS value) in Fig 7D shows that the largest maxF increment (up to 1.7 times) is achieved for the forces of two types of MUs (FR, FF) after synchronization with Method 4 or Method 1, Δt = ± 6 ms. Considering the whole muscle, the observed maximal increment of maxF was about 1.5 times using Method 4, Δt = ± 6 ms. This relative change of maxF was found to be smaller than the force amplitude variances estimated above by the other two force parameters (rangeF and rmsF). This could be explained by the fused nature of

PLOS COMPUTATIONAL BIOLOGY
Synchronization of motor units in a model of the rat medial gastrocnemius muscle maximum force, accumulating both the force mean and variance, the former shown above to be unaffected by synchronization.
• VAF: The synchronization had an important effect on the force variance estimated by VAF, which was decreasing from 100% (NS), greater than 80% (weak synchronization with Δt = ± 2 ms) to down below -200% (strong synchronization with Methods 1, 2 and 4, Δt = ± 6 ms), as shown in Table 3 and Fig 7E. Negative values of VAF are not typical in literature, although Okada et al (2017) [40] reported that VAF decreased from positive to negative values when the analysis time window was widening. The negative VAF in our case was computed due to high variance of the force differences before and after synchronization, which was more than twice larger than the reference variance of the non-synchronized forces. This undoubtedly indicated that the synchronized forces had a very large variance, and they were completely non-synchronized with the forces of the non-synchronized model. Specifically, the synchronization of two types of MUs (S, FR) was noted as most influential on force VAF decrement ( Table 3, Fig 7E), which corresponds to the above observations for largest rmsF values for the output forces of those types of MUs (Fig 7B).
• Force mean spectral frequency: In this context, the synchronization had an important effectdecreasing the meanfreq value by more than 10 Hz, which drops it from 35.6 Hz (NS) down to 24.4 Hz for the muscle force after synchronization with Method 4, Δt = ± 6 ms ( Table 3) It is worth to stress that unlike S or FR MUs, within FF MUs greater synchronization effects were revealed for the weaker subpopulation of units (FF32-FF44) than for the stronger ones (FF45-FF57) in three applied simulations methods (Method 1, 2, 4, Fig 5).

Correlation of the force variance and MU synchronization
The results presented in this section aim to answer the general question of whether the provided synchronization methods regularized by widening the time window for synchronization (Δt = ± 2, ± 4, and ± 6 ms) led to consistent increases in both the level of MU synchronization (CISI) and the induced changes in the force output. Thus, the force parameters, which were most closely correlated to the synchronization design in Methods 1 through 4, could be deduced. The results in Table 4 establish the correlations between the curves in Fig 7A for the level of MU synchronization in the function of Δt (CISI = f(Δt)) and each of the curves in Fig  7B, 7C, 7D and 7F for the trends of the six force parameters as a function of Δt (meanF, rmsF, rangeF, maxF, VAF, meanfreq). The correlations were estimated in the range [−1;+1], where +1 and −1 stand for strongly correlated curves that were directly or inversely proportional, respectively. The results show strong correlations of all parameters, which account for increasing the force variance while increasing the level of synchronization (i.e. rmsF, rangeF, maxF and VAF with an average correlation of about ± 0.97 in all types of MUs). The force mean spectral frequency was indeed inversely proportional to the synchronization level, with an average correlation coefficient of −0.89. The parameter related to the mean force (meanF) was the one least dependent on the synchronization, with an average correlation coefficient of 0.53.

Discussion
There are two different approaches one could use to investigate the synchrony between different MUs and its influence on the developed muscle force. The first one involves assessing experimental recordings of electromyographic signals using needle or surface electrodes and decomposing these signals into individual action potentials [4,[41][42][43][44]. However, the disadvantages of this approach include that only a portion of the active MUs is recorded, it is not possible to distinguish fast from slow MUs and the measured muscle force reflects the force of all active MUs, and even MUs of other muscles. The second method is based on pure modeling, wherein models of the muscle are composed using different MUs [22,32]. These models are based on the Fuglevand et al. approach [31] which contain 120 MUs. However, these authors did not divide MUs into different types (S, FR, and FF). Moreover, the function used for describing the twitch was based only on two parameters: force amplitude and contraction time. The model used in the current paper, is constructed based on experimental data concerning MU twitch and tetanus properties as well as motoneuronal excitability, and has been fully described previously [34]. Here, the experimentally measured twitches are modeled by a sixparameter function and the summation of the twitches into tetanus is established by an experimentally verified mathematical algorithm. In the adopted basic model, it was proven that the firing of all MUs is asynchronous. Then, synchronization was imposed in this basic MU firing arrangement, changing the pattern of pulses of MUs during the steady state of the muscle force using several simulated situations (i.e., four modes of synchronization with the three time windows ± 2, ± 4, and ± 6 ms). In this way, broad investigation of the influence of the synchrony of the three types of the MUs on the developed muscle force and cumulative forces of MUs from the three groups could be performed. The results based on the two used coefficients corMU and CISI showed that the range, the maximum, and the root-mean-square of the forces rose with increased synchronization, while the mean forces remained nearly unchanged. This increase was stronger for fast MUs; notably, these units are mostly responsible for the force instability (muscle tremor) in the context of moderate or strong muscle contractions, wherein fast MUs are recruited into activity.

Models of MU synchronization
To increase the degree of synchronization and to analyze its effects on the muscle and MU forces, we considered the synchronization of pulses of pairs of MUs in the time windows ± 2, ± 4, and ± 6 ms. It is known that synchronization is an effect of a common excitatory input to several motoneurons and that synchronic excitatory postsynaptic potentials (EPSPs) evoked in several motoneurons increase the probability of the simultaneous occurrence of their action potentials [45]. The size of the time windows is related to the duration of EPSPs in rat Table 4. Correlation coefficients between CISI and the six force parameters meanF, rmsF, rangeF, maxF, VAF and meanfreq. The strength of the correlation is coded with a color gradient, highlighting the strong positive (> 0.8) (dense red) and strong negative (< −0.8) (dense blue) correlations.  [46]). Additionally, the applied method resulted in a narrow peak in the cross-interval histogram (Fig 6), similar to that reported for human muscles by De Luca et al. [4], as is typical for short-term synchronization (i.e., the peak centered about zero-time delay 0.5 ± 2.9 ms) and with an average width of 4.5 ± 2.5 ms. For all four proposed modes of synchronization used in the present study, the same range of time windows was applied. The largest (± 6 ms) time window increased the CISI by about 1.5 times for Method 1, about 2.5 times for Methods 2 and 3, and more than three times for Method 4 (see Table 2). The range of differences in the obtained synchronization is similar to that of differences in the CISI reported for trained and nontrained subjects (more than two times higher in weightlifters), changes resulting from conditioning exercise (about 2.5 times higher after the exercise), and those between dominant and nondominant hands (1.6 times higher in the no dominant hand) [43]. The proposed method of inducing synchronization within time windows Δt of variable duration appeared to be an efficient tool in the four tested simulations. For all four methods of synchronization, values of the investigated parameters, which account for the force variance (rmsF, rangeF, maxF and VAF) which characterized amplified force oscillations, changed along with increases in the time window Δt, i.e., when the synchronization degree was augmented (Fig 7). Notably, this change appeared strongest with Method 4 and weakest with Method 3. Meanwhile, the highest value of corMU (74.6) was achieved for Δt = ± 6 ms for FR MUs (Table 2). Moreover, except for in Method 3, the highest values of corMU were observed for FR MUs (Table 2). This observation is surprising in light of previous physiological experiments concerning force decreases/increases resulting from the prolongation/shortening of one IPI during the unfused tetanic contraction ascertained using MUs of the rat medial gastrocnemius muscle [29]. Namely, relative force fluctuations noted for FF and FR MUs were similar and one could expect no differences to exist between these two types of fast MUs in the present simulation study. This methodological approach resulted in the highest synchronization for Method 4 and is reflected by the parameters corMU and CISI in Table 2. It should be stressed that the four methods led to similar effects on muscle force-that is, greater maximal force and higher fluctuations around a mean force-and these increases concerned all three types of MUs, although it should be stressed that this result was obtained for the maximum excitation signal, i.e., a simulation of a very strong contraction, when all MUs were active.

Force
The induced synchronization patterns have random character. First of all, impulses of each individual MU, firing with a given mean frequency, were distributed randomly by an algorithm changing successive IPIs. Second, synchronization was induced only when impulses of two MUs occurred in a close time interval (of 2, of 4 or of 6 ms), which was also a random situation since the time with which the impulse is shifted is random. The synchronization procedure modified the distribution of IPIs for all MUs included in the model and shapes of histograms (S1 Fig in the supplementary material) for non-synchronized and synchronized IPIs (Methods 1-4) were similar to those reported by authors studying MU firing rates in human muscles during voluntary activity [47,48]. Nevertheless, the increased minimum to maximum range of IPIs observed in the histograms in result of four applied methods of synchronization turns our attention to a significance of rate coding. In a recent study on insects Putney et al. (2019) [49] reveled that timing of impulses between motoneuron discharge patterns plays a higher role than a number of pulses, but it is well known that changes in IPIs are crucial for the force development also in the mammalian (cat and rat) muscles [29,30,35,50,51,52].
We decided to induce synchronization within each type of MUs (S, FR or FF, but not randomly between all types of MUs), which corresponded to the primary goal of the study concerning the role of synchronization of particular types of MUs. Moreover, assumptions of the common-drive hypothesis regarding common input to numerous motoneurons in a pool innervating a muscle have some limitations, and there may be several basic patterns of synaptic input organization to motoneurones within a given MU pool. Both peripheral and descending inputs to slow and fast motoneurons differ with respect to latencies or amplitudes of postsynaptic potentials. For example, monosynaptic EPSPs from muscle spindles are the largest in type S MUs, smaller in type FR, and smallest in type FF units [53], disynaptic IPSPs facilitated by rubrospinal volleys are the highest in slow and lowest in fast MUs [54], slow and fast forelimb motoneurons receive different inputs from propriospinal neurons conveying pyramidal commands [55]. All these observations indicate that groups of motoneurons of MUs of each type form specific subpopulations within a motor nucleus.

Effects of synchronization on MU and muscle forces
The influence of the increasing synchronization level on the mean as well as on the maximum force of particular MU types and of the whole muscle was, in general, very weak (i.e., the maximum force increased by up to 5% for the whole muscle and up to 7% for FF MUs), regardless of the synchronization method applied in the model. This confirms the results of previous studies, which also demonstrated that the magnitude of force output and the average force of the muscle were not altered considerably due to synchronization [22,43]. However, an increase in the synchronization time window from ± 2 to ± 6 ms in all cases correlated with a rise in the force of each MU type, with the change being the greatest for synchronization using Methods 1 and 4. Moreover, the present study has revealed certain differences between MU types. Not only did the absolute force increase but also the relative force increased after synchronization; further, they were always the highest for fast MUs (FF and FR) and the lowest for slow MUs. This also confirms previous observations that synchronization may be beneficial during the performance of contractions where rapid force development is required, for which fast MUs should be recruited [18].
On the other hand, it was already mentioned that a muscle can produce smooth contractions due to asynchronous discharges of motor neurons [18,24] and that synchronization increases the variability in the muscle force [22]. Indeed, simulated contractions in our model have confirmed that synchronization substantially influences the range of force oscillations during the steady state of the muscle contraction and the min-max range of modeled forces gradually rose with the increase in the time window for synchronization in each method. This can be partly explained by previous computer simulations indicating that synchronization leads to an increase in the estimated twitch force and to a decrease in the estimated contraction time of an MU [27]. Obviously, absolute values of the min-max range of the force were the lowest for the weakest S MUs, but the ratio of parameters accounting for the force variance (rmsF, rangeF, VAF) between synchronized and NS models was the highest for S MUs for all methods-except Method 4, in which MUs of the same type were synchronized according to the first MU in the group (see Fig 7B, 7C and 7F).
A 100% excitation signal (corresponding to a very strong muscle contraction) used in this model was applied to ensure activation of all MU types, which helped us to elucidate the contributions of the three types of MUs to muscle tremor, which are dependent on the force level [56]. According to the size principle, at a lower excitation signal, a contribution of high-threshold fast MUs (especially those of the FF type) to the force development would be smaller or recruitment would be restricted to low-threshold (S or FR) MUs. The lowest relative force oscillations were noted in FF MUs for all methods of synchronization (Fig 7B and 7E). This observation indicates that slow MUs have the strongest and FF MUs have the weakest relative influence, respectively, on force fluctuations described as muscle tremor and thus partly explains why tremor is best visible during weak contractions, when predominantly slow MUs are recruited. Surprisingly, the model indicated that weaker and stronger FF MUs revealed different synchronization effects (Fig 5). The model is based on data concerning contractile properties of a set of MUs studied in electrophysiological experiments, selected from a large population of units to form a group of MUs reflecting a real number and proportion of MU types (S, FR and FF) in the medial gastrocnemius muscle. Moreover, they were selected to represent mean values of their contractile properties as close as possible to those observed in the whole population. MUs are recruited in the model in order of the increase in the input to motoneurons, and according to the size principle, i.e. from the weakest to the strongest MU within each type. As indicated in Table 1, FF MUs are ordered according to their increasing twitch force and a half of weaker ones (FF32-FF44) have somewhat different properties than the stronger ones (FF45-FF57). Namely, the weaker FF MUs have shorter twitch time parameters than the stronger FF units (the contraction time 13.4±1.7 vs 14.2±1.4 ms, the half-relaxation time 25.4 ±4.6 vs 27.8±4.5 ms), while their firing rates predicted by the model are higher (57.9±14.1 vs 48.0±9.6 ms, respectively), because the firing rate is dependent on the twitch contraction time [34]. Lower firing rate decreases a number of action potentials in the studied time window and most probably provides an explanation of the observed sub-grouping observed in Fig 5. Nevertheless, a question whether stronger FF MUs (those with the highest recruitment threshold) reveal weaker synchronization is open and needs to be studied in separate physiological experiments.

The influence of synchronization on the spectral frequency of the muscle force
To our knowledge, the parameter meanfreq of the force has not been analyzed in muscle modeling in connection with the synchronization of MU firing to date. It should be noted, however, that the power spectral analysis of tremor in the first dorsal interosseous muscle revealed three frequency peaks occurring at around 10 Hz, 20 Hz, and 40 Hz [25], which correspond to our findings concerning the mean spectral frequencies of S, FR, and FF MUs, respectively (Table 3). McAuley et al. [25] concluded that their results reflected the synchronization of MUs at frequencies determined by oscillations within the central nervous system; however, our findings suggest that the force oscillations related to three types of MUs likely contribute to those frequency peaks.
A decrease in the meanfreq was observed in parallel with an increase in the degree of synchronization in all four applied methods. It should be stressed that the mean firing frequencies of all MUs remained unchanged during simulations, due to a constant number of pulses in the analyzed time window (2000 ms). A decrease in force spectral frequencies paired with the occurrence of slower force oscillations. This observation at increased synchronization levels indicates that the force-frequency spectrum depends upon the temporal distribution rather than on the mean firing frequencies of MUs and this conclusion concerns all three types of MUs, despite considerable differences in the meanfreq between S, FR, and FF MUs. The decrease in the meanfreq could not be linked to muscle fatigue, which was not modeled, and should instead be connected with processes of summation of twitches into tetanic contractions.
McAuley and Marsden [57] in their review argued that the physiological tremor in humans is likely of multifactorial origin, with contributions from the 10-Hz range of oscillatory activity of the central nervous system, MU discharge frequencies, reflex loop resonances, and mechanical resonances. However, it must be emphasized that the present results were obtained using the model of a rat muscle, so it is risky to directly compare the frequencies related to different types of MUs collected herein to human data, most of all because rat MUs demonstrate considerably faster contractions and have higher discharge frequencies.

Significance of the results
This is the first study on functional consequences of increasing synchronization in MUs' firing related to the MUs' type. It should be stressed that data concerning the MUs' synchronization have been obtained so far predominantly in human experiments from EMG recordings from a voluntary active muscles, subsequently decomposed into trains of action potentials of MUs. It was found that the synchronization varied across different skeletal muscles [3,15], and increased as an effect of training [18,20,22] or fatigue [23]. It is well known that a certain proportion of different types of MUs is a characteristic feature of a given muscle (which can be composed predominantly of fast or slow MUs) [58,59], which is not permanent during altered motor activity, e.g. the endurance training can evoke adaptation in contractile properties mainly in FR MUs [60], whereas weight-lifting training can modify properties of FF MUs, as well [61]. On the other hand, development of fatigue is predominantly an effect of activity of FF MUs which have low fatigue index [59]. Therefore, the present observations made for separate MU types contribute to better understanding of consequences of MU synchronization in muscles which are different with respect to the MU content and activity levels.

Limitations of the study
The synchronization was studied within a limited interval of 2000 ms, including limited MU stimuli. Their average number was 100 (S MUs), 111 (FF MUs) and 140 (FR MUs) defining an average influence of the synchronization of a single pulse in the range from 0.7% to 1%. Although our analysis of synchronization did not account on this value, it might indirectly determine the resolution for computation of the correlation coefficient and CISI. More detailed analysis of intervals would improve the resolution; however, it is unlikely to change the observed global results, which revealed differences in presence and absence of synchronization much above 1% (Table 2).
It is assumed in the applied model that forces of all co-active MUs reveal algebraic, linear summation, and this is a simplification for several reasons. First, we know that in the modeled muscle (medial gastrocnemius) processes of summation of MUs forces are not linear, and most probably dependent on interactions between co-active MUs, and effects of summation of MU forces progressively decrease with the increased number of active units [62]. Second, the force production in a muscle depends on the muscle length (stretch) [63]. In the present paper we have used for modeling twitch properties obtained under isometric conditions, and therefore the presented results should be considered as modeling of the force produced during contractions of a muscle kept at a constant length. Third, in the present model MUs are not positioned within the muscle, and one should be aware that the interrelationships between muscle fibers of co-active MUs influence the cumulative force, and overlapping of MU territories most likely negatively influence the muscle force output [62]. A problem of interrelationships between neighboring MUs deserves to be modeled, although till now we have a limited input data to propose a reliable model.

Conclusion
The present study revealed that, regardless of the method used for the synchronization of MU firings, the increase in the synchronization index had a negligible effect on the mean force of the developed contractions yet influenced muscle tremor by increasing force oscillations and further highlighted that these results were observed for all three types of MUs. A parallel decrease in the mean spectral frequency of the force indicated that, in the synchronized models, the force oscillations were slower despite higher magnitudes. The synchronization of fast MUs led to higher increases in the range of the force variability and the force root-meansquare in comparison with that of slow MUs. On the other hand, relative changes in the latter parameters in the synchronized simulations were the highest for slow MUs, indicating their significant contribution to muscle tremor, especially during weak contractions.