SYDE544 -Assignment #4 - Solved

setting up the stage 
This part of the work, we need to setup the stage, i.e. generate a number of parameters pertaining to the characteristics of a pool of motor units (MU) that will be activated in a reasonably physiological manner. Then we need to put them to work, i.e. generate the timing of the firings of these MUs.  

In this simulation, we consider the investigation of a motor unit pool with an unknown number of MUs. However, for surface recordings, we don’t usually have access to the activities of all the Mus in the pool. Rather, we assume we have access superficial 120 MUs. Since there is no spatial organization order with respect to the recruitment thresholds of these units, it is reasonable to assume these 120 MUs are representative to the physiological recruitment order: smaller units are recruited first (at lower threshold), larger units are recruited at higher threshold. Note that the size here refers to the size of the axon of the motoneuron and the size of force twitch generated by the MU, not the size of the action potential detected at skin surface. This is because a deep unit (very small action potential at surface) could be a large unit with large force output. Conversely, a superficial unit (very large action potential at surface) could be a smaller unit with small force output.

Let’s define the recruitment threshold of the ith MU as:

 where the constant a is:
  
(1)
  (2)

 

where RR is the largest recruitment threshold of MU pool (the recruitment of the largest unit). This means once the neural drive applying to the MU pool is larger than  , then the ith MU will be recruitment, and start to firing. With this simplified, yet rather realistic modelling, we can establish the ordered recruitment sequence of the MU pool: a majority of the MUs have small recruitment thresholds, and a small number of MUs have a high threshold. In reality, it is physiologically reasonable to have a 30-fold difference between the highest threshold and the smallest threshold. Therefore, we can set RR=30 (arbitrary unit) in Eqn. (2).

Once the ith MU is recruited at its recruitment threshold, RTEi, it would fire at its minimal firing rate. And, as the neuronal drives increases, it would increase their firing rates. A simple but realistic model for this behavior (rate coding) is:

                                                                                                (3)

where   is the firing rate of the ith MU in response to the applied neural drive to the MU pool E, which is common for all MUs. The coefficient   is a gain factor for the ith MU, and 8 (Hz) is the minimal firing rate (at the recruitment threshold), which is assumed to be identical for all units. To further simplify our simulation, let’s further assume the maximal firing rate of all units to be 35 Hz, and the maximal neural drive is 36 (1.2 times larger than the maximal recruitment threshold). As such, the gain factor can be found by:

                                                                                                           (5)

Problem 1a: Obtain the   for each of the 120 MUs. (complete the code at Line 15-21 of the script.). (0.5 pt)

Problem 1b: Obtain the   for each of the 120 MUs (complete the code at Line 24-30). (0.5 pt)

 

Part two: Simulating EMG 
 Problem 2, Using  and  in Problem 1, generate a realization of the firing timings of all 120 MUs over approximately 10s, at each of the neural drive levels:

  

When generate the firing timings, the mean firing rate   should be obtained from Eqn (3), the inter-firing interval should be modelled as a Gaussian variable with mean of   and coefficient variation of 15%. (Complete the code at Line 33-39). (2 pt)

Problem 3: In the provided data file, you can find a variable named, MUAPs. It contains the MUAPs of the 120 MUs (each column for a MUAP). The sampling rate is 4096 sample per second. If you plot these MUAPs, you should see that there is no relationship between the size of the MUAP and MU id, as explained earlier. Using the firing timings generated in Problem 2, simulate the MUAP trains of all the units at the 9 levels of neural drive. Add the MUAP trains together, you would obtain simulated EMG signals at these 9 levels of neural drive. (Complete the code at Line 42-50.) (1 pts)

 

Problem 4: Calculate the SNR of EMG, according to the formula discussed during our lecture. You can take a 2-s portion of the middle of the 10-s long simulated data, avoiding the beginning and the end parts. (Complete the code in Line 52-55) (2 pt)

 

Part three: Simulating Force output 
Problem 5: In the provided data file, you can also find another variable named, F. It contains the force twitches of the 120 MUs (each column for a MUAP). For each MU, a firing will generate one force twitch. The sampling rate is also 4096 sample per second. If you plot these force twitches, you should see that the force twitches are nice ordered: the twitch of the first unit is the smallest and twitch of the last unit is the largest. Using the firing timings generated in Problem 2, simulate the force output of all the 120 MUs, at the 9 levels of neural drive. The summation of the force output of all unit is the overall force output. (Complete the code in Line 58-62) (2 pts)

Problem 6: Calculate the SNR of force output, with a definition similar to the one for EMG, i.e. the square of the mean force output divided by the variance of the force output.

Ideally, you should use the same time period as in the SNR calculations of EMG.

(Complete the code in Line 64-66). (1 pt)

 

Problem 7: plot the SNR of both EMG and force output against the level of neural drive, and comment on your observation. (1 pts)

 

Reading material  
 

A more comprehensive simulation methods, on which the current assignment is based, can be found in the following paper (attached):

 

AJ Fuglevand, DD Winter, and AE Patla, “Models of Recruitment and Rate Coding

Organization in Motor-Unit Pools”, Journal of Neurophysiology, Vol. 70, No. 6, pp. 24702487, Dec. 1993.