07/11

b_factory.py

@@ -10,12 +10,42 @@
 import scipy as sp
 
 # import the necessary functions
-from s_functions import simulate_neuron, simulate_nm_conc
+from s_functions import simulate_neuron, simulate_nm_conc, simulate_fluorescence_signal
 
 
 # define the background fluorescence due to the tissue
 bg_tissue = 1.5
 
+
+# functino to generate heatmap
+def bleach_dnm_heat(tau_values, nm_conc_input, var):
+
+    # create an array to store the snr values
+    snr = np.zeros((tau_values.size, tau_values.size))
+
+    # BONUS: find way to do it that's more effecient -- this is O(n^2)
+    for i in range(len(tau_values)):
+        for j in range(len(tau_values)):
+            progression, progression_sub = simulate_fluorescence_signal(nm_conc=nm_conc_input,variance=var, tau_d=tau_values[i],tau_nm=tau_values[i], tau_tissue=tau_values[j])
+            snr[i,j] = np.abs(np.mean(progression_sub[3])/np.std(progression_sub[3]))
+            #print(np.std(signal[3])) -- for testing
+
+    # use log tau values to make a simple scale
+    log_tau_values = np.log10(tau_values)
+
+
+
+    # Generate the heatmap
+    start, stop = log_tau_values[0], log_tau_values[-1]
+    plt.imshow(snr, cmap='magma', extent=[start, stop, stop, start])
+    plt.colorbar(label='snr')
+    # plt.xlabel('bleach time constant in tissue fluorescence')
+    # plt.ylabel('bleach time constant in dye + nm flourescence')
+    plt.title('variance = {}'.format(var))
+    #plt.show()
+
+    return snr
+
 # ANALYSIS 1: bleaching factor acting on the [NM] contribution to F
 def bleach_nm(K_D, tau, F_max, F_min, nm_conc, bline_len=5000):
 
@@ -120,37 +150,37 @@ def bleach_dnm(K_D, tau, F_max, F_min, nm_conc, bline_len =5000):
         delta_f[i] = f_sub[i]-f0
         delta_ft_f0[i] = delta_f[i]/(f0)
 
-    print(f0)
+
 
     # plot f, f - fit, df then df/f
-    plt.figure()
-    plt.subplot(2,2,1)
-    plt.plot(t,f, label='f')
-    plt.xlabel('time (ms)')
-    plt.ylabel('f')
-    plt.legend()
-
-    plt.subplot(2,2,2)
-    plt.plot(f_sub, label='f subtracted')
-    plt.plot(fit, label='fit')
-    plt.xlabel('time (ms)')
-    plt.ylabel('f-exp')
-
-    plt.legend()
-
-    plt.subplot(2,2,3)
-    plt.plot(t,delta_f, label='df')
-    plt.xlabel('time (ms)')
-    plt.ylabel('delta f')
-    plt.legend()
-
-
-    plt.subplot(2,2,4)
-    plt.plot(t,delta_ft_f0, label = 'df/f')
-    plt.xlabel('time(ms)')
-    plt.ylabel('Delta F/F0')
-    plt.tight_layout()
-    plt.legend()
+    # plt.figure()
+    # plt.subplot(2,2,1)
+    # plt.plot(t,f, label='f')
+    # plt.xlabel('time (ms)')
+    # plt.ylabel('f')
+    # plt.legend()
+
+    # plt.subplot(2,2,2)
+    # plt.plot(f_sub, label='f subtracted')
+    # plt.plot(fit, label='fit')
+    # plt.xlabel('time (ms)')
+    # plt.ylabel('f-exp')
+
+    # plt.legend()
+
+    # plt.subplot(2,2,3)
+    # plt.plot(t,delta_f, label='df')
+    # plt.xlabel('time (ms)')
+    # plt.ylabel('delta f')
+    # plt.legend()
+
+
+    # plt.subplot(2,2,4)
+    # plt.plot(t,delta_ft_f0, label = 'df/f')
+    # plt.xlabel('time(ms)')
+    # plt.ylabel('Delta F/F0')
+    # plt.tight_layout()
+    # plt.legend()
 
 
     return delta_ft_f0

f_rate2.py

@@ -0,0 +1,140 @@
+# This script is going to be used to compare the df/f vs firing rate plot for different bleaches
+
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+import pandas as pd
+import scipy as sp
+
+from s_functions import simulate_neuron, simulate_nm_conc, simulate_fluorescence_signal
+from b_factory import bleach_dnm_heat
+
+
+# define the base function
+
+def dff_v_activity(neuron_activity, progression_array, window=100):
+
+    number_blocks = int(progression_array[3].size/window)
+    signal_points = []
+
+
+    # array of neural activity
+    firing_rate = []
+
+    for block in range(number_blocks):
+
+        # define the interval
+        start = (block*window)
+        stop = ((block+1)*window+1)
+
+        # get average signal for 100 ms sliding window
+        signal_average = np.average(progression_array[3][start:stop])
+
+        # store it
+        signal_points.append(signal_average)
+
+        # Check the activity of the firing neuron for current window
+        spikes = np.size(np.nonzero(neuron_activity[start:stop]))
+
+        # convert it to Hz and add it to the firing rate array
+        firing_rate.append(spikes*1000/100)
+
+
+    # get the unique firing rates
+    unique_frates = np.unique(firing_rate)
+
+    # array to store the average dff values
+    average_signal_points = []
+
+    # convert the list to an array
+    signal_points = np.array(signal_points)
+
+
+    # average the df/f values at each of the unique values
+    for i in range(len(unique_frates)):
+
+        # get the indices where the firing rate was a given value, x
+        indices = np.where(firing_rate==unique_frates[i])[0]
+
+        # average all the corresponding df/f values
+        ave_dff = np.average(signal_points[indices])
+        average_signal_points.append(ave_dff)
+
+
+    # fit the average signal points 
+    poly = np.polyfit(unique_frates,average_signal_points,1)
+    fit = np.polyval(poly,unique_frates)
+
+    return unique_frates, average_signal_points, fit
+
+
+def plot_dff_v_activity(firing_rates, dff, fit):
+
+
+    plt.plot(firing_rates,dff,'o')
+    plt.plot(firing_rates,fit, label='fit')
+    plt.xlabel('Firing rates(Hz)')
+    plt.ylabel('Average df/f signal')
+    plt.title('Signal vs activity plot (average)')
+    plt.legend()
+    plt.show()
+
+
+# Test 1
+
+# At different bleaches simulate a neuron firing and get the signal vs firing rate plot
+# bleaches = np.logspace(5,7,5)
+
+# # generate a firing neuron
+# neuron = simulate_neuron(n_timesteps=70000,firing_rate=13)
+
+# # generate nm_conc 
+# nm_conc_input, nm_b_conc, nm_r_conc = simulate_nm_conc(neuron,nm_conc0=0,k_b=0.6, k_r=0.4,gamma=0.004)
+
+# # plot the df/f vs firing rate at diff bleaches
+# plt.figure()
+
+# for i in range(len(bleaches)):
+
+#     # genrate the progression array for fluorescence
+#     progression, progression_sub = simulate_fluorescence_signal(nm_conc=nm_conc_input,tau_d=bleaches[i],tau_nm=bleaches[i], tau_tissue=10e9)
+
+#     # plot the signal vs firing rate
+#     x, y, fit = dff_v_activity(neuron,progression_sub)
+#     #plt.scatter(x,y,label=np.round(bleaches[i]))
+#     plt.plot(x,fit,label=np.round(bleaches[i]))
+
+
+# plt.xlabel('Firing rates(Hz)')
+# plt.ylabel('Average df/f signal')
+# plt.title('Signal vs activity plot')
+# plt.legend(title='bleach time constants')
+# plt.show()
+
+
+# Test 2:
+
+# check the effect of changing the variance of the gaussian noise on ftissue on the snr 
+var_values = np.array([0.0001,0.001,0.01,0.1,1,3])
+var_v = np.array([1,3,10])
+
+# bleach time constants for heatmap
+specific_taus = np.logspace(5,7,20)
+
+# generate a firing neuron
+neuron = simulate_neuron(n_timesteps=70000,firing_rate=13)
+
+# generate nm_conc 
+nm_conc, nm_b_conc, nm_r_conc = simulate_nm_conc(neuron,nm_conc0=0,k_b=0.6, k_r=0.4,gamma=0.004)
+
+# for different variances, get the heatmap
+plt.figure()
+for i in range(len(var_v)):
+    plt.subplot(3,2,i+1)
+    bleach_dnm_heat(specific_taus,nm_conc_input=nm_conc, var = var_v[i])
+
+plt.suptitle('SNR vs bleach strength at different variance for ftissue', size = 16)
+plt.tight_layout()
+plt.show()
+
+

s_functions.py

@@ -29,21 +29,13 @@ def simulate_neuron(n_timesteps, firing_rate, number=1):
     firing_neuron = firing_neuron.astype(int)
 
 
-    # then make a plot of it!
-    # plt.plot(firing_neuron)
-    # plt.xlabel('timesteps')
-    # plt.ylabel('Neuron activity')
-    # plt.title('Neuron Activity ({}Hz) over {} timesteps'.format(firing_rate,n_timesteps))
-    # plt.show()
-
 
     # check exactly how many spikes were produced: to see if it works
     n_spikes = np.size(np.nonzero(firing_neuron))
 
     # print simulated neuron summary:
     #print('Simulated neuron with {} spikes in {} timesteps ({} Hz).'.format(n_spikes, n_timesteps, firing_rate))
 
-
     return firing_neuron
 
 # test 1
@@ -153,10 +145,10 @@ def plot_nm_conc(nm,start,stop,colour='b', plotlabel = ''):
 
 
 
-def simulate_fluorescence_signal(tau_d, tau_nm, tau_tissue, nm_conc, K_D = 1000, F_max = 45, F_min = 10, bline_len=5000):
+def simulate_fluorescence_signal(tau_d, tau_nm, tau_tissue, nm_conc, variance=0.0005, K_D = 1000, F_max = 45, F_min = 10, bline_len=5000):
 
-    # autofluorescence
-    f_tissue = np.random.normal(0.002,0.0005,nm_conc.size)
+     # noisy autofluorescence
+    f_tissue = np.random.normal(0.002,variance,nm_conc.size)
 
     # create timesteps 
     n_timesteps = nm_conc.size
@@ -170,16 +162,8 @@ def simulate_fluorescence_signal(tau_d, tau_nm, tau_tissue, nm_conc, K_D = 1000,
     # calculate F: derived from eq 2 in Neher/Augustine
     f = bleach_tissue*f_tissue + (bleach_d*K_D*F_min + bleach_nm*nm_conc*F_max)/(K_D + nm_conc)
 
-    # calculate f0 by getting the median value of the bottom 70% of previous f values
-    percentile_mark = np.percentile(f,70)
-    f0 = np.median(f[f<percentile_mark])
-
-    # df calc -- median value method
-    df = f-f0
-    df_f_med = df/f0
 
-
-    # SUBTRACTED VERSION: fit an exponential to remove the bleaching trend 
+    # fit an exponential to remove the bleaching trend 
 
     # define an exponential function that we'll use as the basis for the fit
     def exp_decay(t,a,b):
@@ -193,13 +177,23 @@ def exp_decay(t,a,b):
 
     # define the fitted function
     fit = exp_decay(t,a_fit,b_fit)
-
+    
     # subtracted f
     f_subtracted = f - fit
 
+
     # to correct for negative values in the fluorescence that result in a -ve df/f
     f_alt = f_subtracted + np.max(np.abs(f_subtracted))
 
+
+    # calculate f0 by getting the median value of the bottom 70% of previous f values
+    percentile_mark = np.percentile(f,70)
+    f0 = np.median(f[f<percentile_mark])
+
+    # df calc -- median value method
+    df = f-f0
+    df_f_med = df/f0
+
 
     # df/f with the subtracted formula
     # subtracted signal
@@ -233,11 +227,12 @@ def exp_decay(t,a,b):
 
         # average value
         df_f_ave[i] = (f[i] - f0_ave)/f0_ave
-        f0_averages[i] = f0_ave
+        f0_averages[i]=f0_ave
 
     # define progression arrays for f, df, df/f
     progression = []
     progression_sub = []
+    #progression_multi = []
 
     # without the subtracted exponent
     progression.append(f)