In a recent preprint (McKiernan, 2015), I describe the results of experiments to test the effects of genetically manipulating motor neuron excitability on locomotor activity in Drosophila larvae. I expressed the recombinant line known as 'Electrical Knock-In' (EKI) (Hartwig et al., 2008) in select motor neurons. EKI expression in motor neurons reduces potassium currents and causes cells to fire soooner, faster, and at lower current levels. In the current study, I used the FLP/FRT system (Golic & Lindquist, 1989) to express EKI in select motor neurons while other motor neurons in the same animal remained wildtype. In this way, both control and experimental conditions could be compared in the same animal. Because recording directly from the motor neurons during spontaneous locomotor activity presents several challenges, and since potentials recorded from target muscles are one-to-one with potentials arising from innervating motor neurons (Choi et al., 2004), I recorded from the muscles as a proxy for neural activity. In each animal, the bursting activity from a muscle innervated by wildtype motor neurons and a muscle innervated by at least one motor neuron (MN1-Ib or MNISN-Is) expressing EKI was recorded, quantified, and compared. For more details on study design and methods, please see McKiernan (2015).
I have made all the data and analysis code associated with this study publicly available in a github repository. In this notebook, I describe how to access the raw muscle recordings, analyze the bursting data, and recreate the tables and figures provided in McKiernan (2015).
To have all plots render within the notebook, enter the following command:
%matplotlib inline
All the raw electrophysiological recordings are available in both ABF (Axon) and SMR (Spike2) file formats. Although both these formats are proprietary, we can access the data and plot the recordings using the python package Neo. Functions for reading and plotting both ABF and SMR files are found in 'readRecordings.py'. Run the file by entering the following command. This will import the necessary python packages and define the functions.
%run readRecordings.py
The neo package appears to rescale the voltage when reading ABF but not SMR files. So for the moment, we will preferably read and plot the SMR files so we can see the raw voltages. Using the 'readSMR' function, we specify the file name, and the minimum and maximum x- and y-axis values for plotting. We can choose to view the entire recording by setting xmin = 0 and xmax to a sufficiently large value (this will depend on the length of the recording). Or, we can choose to plot only that portion of the recording included in the final analysis. Not all electrical activity was included in the analysis due to problems with electrode placement, stability, or a lack of crawling-related rhythmic activity. For example, for the following recording, only bursts between time points ~163 and 276 seconds were analyzed, so xmin and xmax are set accordingly. (To determine the relevant time points for a given recording, please refer to the burst times listed in the master data spreadsheet provided in the repository.) The minimum and maximum voltage values (ymin, ymax) are set at -70 and 0mV, respectively. Thus, to read the file, type the following:
readSMR('SMR-files/09706000.SMR',xmin=163,xmax=276,ymin=-70,ymax=0)
The blue and green recordings are from two different channels, recorded simultaneously from muscles in neighboring abdominal segments. Note that there is variability in both the baseline membrane potential and the amplitude of the voltage signal across channels. This variability was common and due in part to muscle movement during locomotor activity and related differences in stability across recordings. Glass electrodes were pulled to have very long, flexible tips that moved with the muscles during contractions. However, depending on factors such as the angle of electrode entry, some recordings were more stable than others. These factors did not, however, affect the quantification of bursting activity for the recordings analyzed here.
Burst times were extracted from the raw recordings and stored in csv spreadsheets. The master spreadsheet contains information on each animal and the burst times from all recorded channels. Csv files containing subsets of recordings for analysis have also been included for use in the following examples. The functions for data analysis are found in 'burstAnalysis.py'. To run the file, type the following command. This will import packages and define the functions.
%run burstAnalysis.py
Next, we specify the csv files to be analyzed. In each animal, two muscles were recorded in neighboring segments. One muscle was innervated by control motor neurons (wildtype; WT), and the other by at least one motor neuron expressing the genetic manipulation EKI. It is assumed that WT and EKI burst times are stored in separate csv files. Two sample csv files are included: 'MN1-Ib_EKI.csv' contains the burst times for muscles innervated by the motor neuron MN1-Ib expressing EKI and 'MN1-Ib_WT.csv' includes the corresponding WT channels recorded in the same animals. To define the csv files and their respective keys, type:
csvFiles={'EKI':'MN1-Ib_EKI.csv','WT':'MN1-Ib_WT.csv'}
Then, extract the data from the csv files. To do this, we run a short for loop that loops over the keys in csvFiles ('EKI' and 'WT') and extracts the data for each.
data=dict()
for k in csvFiles.keys():
data[k]= getCSVData(csvFiles[k])
In total, there were four animals in this sample. So, the EKI and WT data each consist of four sets of burst times. Note that while the EKI and WT channels from the same animal necessarily have the same number of bursts, burst number varies across animals. To print the burst times for each EKI channel, type:
for n in np.arange(0,len(data['EKI'])):
print data['EKI'][n]
[81.21, 110.02, 112.93, 133.61, 136.51, 161.76, 164.83, 185.34, 187.73, 205.83, 209.97, 229.75, 232.89, 246.32, 250.36, 263.0, 267.46, 283.98, 287.66, 305.07, 309.53, 321.23, 326.37, 338.75, 343.1, 357.1, 362.04, 377.77, 383.23, 398.81, 403.68, 413.12, 417.9, 428.34, 432.9, 444.33, 449.79, 459.65, 469.98, 472.76] [232.17, 236.88, 239.74, 244.72, 248.16, 252.37, 257.24, 261.26, 264.42, 270.29, 273.66, 279.53, 290.34, 294.52, 300.89, 304.98] [326.64, 349.09, 354.93, 368.4, 374.46, 391.07, 397.36, 414.42, 420.14, 433.08, 439.65, 449.98, 455.2, 467.31, 473.05, 482.96, 487.55, 498.72, 504.25, 515.31, 523.76, 532.42, 538.89, 547.76, 554.65, 563.73] [11.67, 18.34, 21.44, 25.91, 28.33, 34.58, 36.95, 43.31, 45.57, 52.24, 54.71, 62.07, 66.85, 73.47, 76.42, 83.62, 86.24, 93.71, 98.02, 104.59, 108.21, 112.52, 118.78, 128.13, 146.95, 151.52, 154.3, 162.34, 171.17, 176.8, 179.58, 187.89, 191.3, 199.92, 204.18, 212.48, 222.26, 226.15, 230.04, 242.02]
To print the burst times for the corresponding WT channels from each animal, type:
for n in np.arange(0,len(data['WT'])):
print data['WT'][n]
[81.38, 109.9, 112.72, 132.9, 136.22, 160.97, 164.37, 184.63, 187.47, 205.04, 209.03, 229.33, 232.63, 246.06, 249.05, 262.85, 267.25, 282.46, 286.61, 304.44, 309.11, 320.86, 324.64, 337.33, 341.58, 356.95, 360.09, 377.3, 382.12, 397.34, 402.42, 412.18, 416.11, 426.97, 431.01, 442.97, 446.7, 458.34, 462.64, 471.45] [231.86, 235.99, 239.35, 244.49, 247.46, 252.6, 256.12, 260.21, 263.88, 269.75, 272.88, 278.87, 289.92, 293.82, 299.96, 304.87] [325.07, 346.4, 354.26, 366.15, 373.11, 388.6, 395.9, 407.91, 418.24, 430.92, 437.43, 448.21, 454.83, 465.83, 472.34, 481.66, 487.05, 497.26, 503.55, 514.21, 523.64, 531.16, 538.23, 546.43, 553.27, 562.36] [11.72, 18.39, 21.44, 25.91, 28.43, 34.58, 37.1, 43.31, 45.72, 52.4, 54.97, 61.86, 67.27, 73.63, 76.42, 83.67, 86.3, 93.97, 97.7, 105.32, 108.16, 113.21, 118.88, 128.34, 146.74, 152.15, 154.25, 162.98, 171.23, 177.64, 179.37, 188.83, 191.15, 200.4, 204.13, 213.06, 222.47, 226.99, 230.25, 242.07]
Now, from the burst times we can calculate bursting measures, such as burst duration, cycle duration, duty cycle, and quiescence interval.
Measures:
Again, since there are four animals in this sample, the EKI and WT data each consist of four sets of burst measures. All measures are expressed in seconds, except for duty cycle which is unitless. To calculate the burst measures, run:
bMeasures=dict()
for k in data.keys():
bMeasures[k]= calcBurstMeasures(data[k])
To illustrate the results, we can print all the burst durations for the EKI channels.
for n in np.arange(0,len(bMeasures['EKI']['burstDur'])):
print bMeasures['EKI']['burstDur'][n]
[ 28.81 20.68 25.25 20.51 18.1 19.78 13.43 12.64 16.52 17.41 11.7 12.38 14. 15.73 15.58 9.44 10.44 11.43 9.86 2.78] [ 4.71 4.98 4.21 4.02 5.87 5.87 4.18 4.09] [ 22.45 13.47 16.61 17.06 12.94 10.33 12.11 9.91 11.17 11.06 8.66 8.87 9.08] [ 6.67 4.47 6.25 6.36 6.67 7.36 6.62 7.2 7.47 6.57 4.31 9.35 4.57 8.04 5.63 8.31 8.62 8.3 3.89 11.98]
We can also print all the burst durations for the corresponding WT channels.
for n in np.arange(0,len(bMeasures['WT']['burstDur'])):
print bMeasures['WT']['burstDur'][n]
[ 28.52 20.18 24.75 20.26 17.57 20.3 13.43 13.8 15.21 17.83 11.75 12.69 15.37 17.21 15.22 9.76 10.86 11.96 11.64 8.81] [ 4.13 5.14 5.14 4.09 5.87 5.99 3.9 4.91] [ 21.33 11.89 15.49 12.01 12.68 10.78 11. 9.32 10.21 10.66 7.52 8.2 9.09] [ 6.67 4.47 6.15 6.21 6.68 6.89 6.36 7.25 7.67 7.62 5.05 9.46 5.41 8.73 6.41 9.46 9.25 8.93 4.52 11.82]
We can do the same for the other three measures by replacing 'burstDur' in the for loops above with 'cycleDur', 'dutyCycle', or 'qI', as necessary.
Now that we have all the measures, we can make a few calculations. First, let's specify a coarse enough partition of the sample space. (Note: This binning is different from that which will later be used for plotting.)
Warning: In a few animals, multiple rhythmic bursting bouts were separated by more than a minute, leading to one or two large cycle durations and quiescence intervals. So if the measure being analyzed is one of these, set 'maxPt' higher (e.g. 200) to avoid returning an error that there are values above the interpolation range.
binSize=1e-3; maxPt= 30.0; binEdges= sc.arange(0,maxPt,binSize)
Next, we calculate the cumulative distribution functions (CDFs) for each measure (burst duration, cycle duration, etc.) for all the animals in the EKI and WT samples.
measureNames= bMeasures['WT'].keys()
measureNames.sort()
cdfs={'EKI': dict(), 'WT':dict() }
for k1 in cdfs.keys():
for k2 in measureNames:
cdfs[k1][k2]=calcCDFs(sampleList=bMeasures[k1][k2],binEdges=binEdges)
Now, average over the CDFS for all EKI or WT channels to get the sample CDFs. The following is an example for the measure burst duration ('burstDur').
meanCDF=dict()
for k in cdfs.keys():
meanCDF[k]= cdfs[k]['burstDur'].mean(0)
Quantiles for the different measures can be found by using Kolmogorov´s theorem and interpolation with respect to the average CDFs. If we extract quantiles by finding the inverse images of 0.25 to 1.00 in steps of 0.25 with respect to the average CDFs, we will obtain the first (q1), second (q2; median), and third (q3) quartiles. To do this, run the following for loops:
getQuants=dict()
for k in meanCDF.keys():
getQuants[k]= sc.interpolate.interp1d(meanCDF[k],binEdges)
quartiles=dict()
for k in getQuants.keys():
quartiles[k]= getQuants[k](sc.arange(0.25,1.00,0.25))
We can also calculate the minimum and maximum quantiles for each sample.
minQuants=dict()
maxQuants=dict()
for k in getQuants.keys():
minQuants[k]= getQuants[k](0.01)
maxQuants[k]= getQuants[k](1.0)
We can then print the minimum, maximum, and quartiles values obtained from the EKI and WT samples. These printouts can be copied directly into a LaTeX file to create tables. This procedure was used to create all tables included in the preprint McKiernan (2015). Note that the printout below lists only the burst durations for each sample. The same procedure is repeated for the other three measures to complete the tables.
str0='min & q1 & q2 & q3 & max \\'
str=dict()
for k in quartiles.keys():
print k
str[k]='%.2f & %.2f & %.2f & %.2f & %.2f \\'%(minQuants[k], quartiles[k][0], quartiles[k][1], quartiles[k][2], maxQuants[k])
print str0,'\n',str[k],'\n'
WT min & q1 & q2 & q3 & max \ 3.90 & 5.87 & 8.93 & 12.01 & 28.52 \ EKI min & q1 & q2 & q3 & max \ 2.78 & 4.98 & 8.66 & 12.64 & 28.81 \
For statistical comparisons between the EKI and WT samples, we can use the quartiles above, or we can extract quantiles over a finer range.
quantiles=dict()
for k in getQuants.keys():
quantiles[k]= getQuants[k](sc.arange(0.01,1.00,0.01))
We can run a Wilcoxon rank-sum test to determine whether the two samples are statistically different.
z,p=sc.stats.ranksums(quantiles['WT'],quantiles['EKI'])
print z,p
0.281530046301 0.778303885984
The test returns a p-value of 0.778, so we accept the null hypothesis that the burst durations from the EKI and WT samples come from the same distribution. We can run the same analysis for the three other measures (cycle duration, etc.).
Now that we have all the data extracted and analyzed, we can visually compare the EKI and WT channels. Let's start by plotting the cumulative distribution functions (CDFs) for each channel recorded in each animal, as well as the group average CDFs.
bins= calcBins()
graphMultiCDF(bins,bMeasures)
The dotted lines represent the CDFs from individual animals, each with an EKI (red) and WT (blue) channel. The solid lines represent the average of all the EKI or WT channel CDFs. It can be seen that the CDFs from individual animals as well as the averaged CDFs largely overlap, indicating no difference between the EKI and WT samples.
For illustration purposes, we can also compare the probability density functions (PDFs) averaged over the EKI and WT channels.
pdfComparison(label1='EKI',label2='WT',maxRF=0.35,alpha1=1,alpha2=1,skinfactor=0.3)
MN1-Ib_WT.csv MN1-Ib_EKI.csv
The results show that the EKI (white) and WT (black) PDFs largely overlap, indicating no difference between the samples on any measure. The figure above is the same as Fig. 3 in the preprint McKiernan (2015). The same procedure described above can be followed to analyze and visualize bursting data from animals in which the motor neuron MNISN-Is expressed EKI to recreate Fig. 5, or for data from animals expressing EKI in both MN1-Ib and MNISN-Is to recreate Fig. 7.