10/11

Bleach_compare.py

@@ -10,8 +10,8 @@
 import scipy as sp
 
 # import functions from the simulation and bleaching libraries
-from s_functions import simulate_neuron, simulate_nm_conc, simulate_flourescence_signal
-from b_factory import bleach_nm, bleach_dnm, bleach_t, bleach_all
+from s_functions import simulate_neuron, simulate_nm_conc, simulate_fluorescence_signal
+from b_factory import bleach_nm, bleach_dnm, bleach_t, bleach_all, bleach_dnm_heat
 
 
 # ANALYSIS 1: 
@@ -35,15 +35,33 @@
 
 
 
+#ANALYSIS 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()
 
-# ANALYSIS 2: comparing how the different bleaches affect the signal to noise ratio
 
-# define the signal to noise ratio
 
-# get the signal t
 
 
 
@@ -54,8 +72,6 @@
 
 
 
-
-
 # PREVIOUS ANALYSIS: comparing the different contributions at one timescale for the bleach factor
 # plt.plot(t,b1,label='bleach 1')
 # plt.plot(t,b2,label='bleach 2')

Bleaching w:out_bline.py → Bleaching G0.py

@@ -1,5 +1,5 @@
-# This script was the original analysis script and
-# has since been updated -- cf b_compare or analysis
+# This script was the original bleaching analysis script and
+# has since been updated -- cf bleach_compare
 
 # it has the bleaching analysis with the df/f computed without the 
 # moving baseline calculation of f0

Simulation Script.py

@@ -1,4 +1,24 @@
-#updates coming -- this will probably be the script we call the functions from
+# This script performs the whole simulation in one go: neuron-->fluorescence signal
 
 
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+import pandas as pd
+import scipy as sp
+from scipy.optimize import curve_fit
+
+# bring the necessary functions
+from s_functions import simulate_neuron, simulate_nm_conc, simulate_fluorescence_signal
+
+# simulate a neuron
+timesteps = 70000
+rate = 13
+firing_neuron = simulate_neuron(n_timesteps=timesteps,firing_rate=rate)
+
+# 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, timesteps, rate))
 

SimulationS.py

@@ -0,0 +1,37 @@
+# This script performs the whole simulation in one go: neuron-->fluorescence signal
+
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+import pandas as pd
+import scipy as sp
+
+# bring the necessary functions
+from s_functions import simulate_neuron, simulate_nm_conc, plot_nm_conc,simulate_fluorescence_signal, plot_f_signal
+
+# FUTURE UPDATE: Take as arguments the timesteps and rate -- from the command line 
+
+# STEP1: simulate a neuron
+timesteps = 70000
+rate = 13
+firing_neuron = simulate_neuron(n_timesteps=timesteps,firing_rate=rate)
+
+# 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, timesteps, rate))
+
+# STEP2: simulate nm dynamics
+nm_conc_in, nm_b_conc_in, nm_r_conc_in = simulate_nm_conc(firing_neuron,nm_conc0=0,k_b=0.6, k_r=0.4,gamma=0.004)
+print('Simulated nm dynamics from neuron')
+plot_nm_conc(nm_conc_in, nm_b_conc_in, nm_r_conc_in)
+
+# FUTURE UPDATE: take in as arguments from the command line the bleach time constant for the dye and nm
+bleach_time = 10e5
+
+# STEP3: simulate fluorescence signal from these dynamics and plot it
+progression_in, progression_sub_in = simulate_fluorescence_signal(tau_d=bleach_time, tau_nm=bleach_time, tau_tissue=10e9, nm_conc=nm_conc_in)
+plot_f_signal(progression_in, progression_sub_in, nm_conc_in)
+print('Simulated fluorescence signal from nm dynamics')
+

s_functions.py

@@ -29,7 +29,6 @@ def simulate_neuron(n_timesteps, firing_rate, number=1):
     firing_neuron = firing_neuron.astype(int)
 
 
-
     # check exactly how many spikes were produced: to see if it works
     n_spikes = np.size(np.nonzero(firing_neuron))
 
@@ -47,6 +46,31 @@ def simulate_neuron(n_timesteps, firing_rate, number=1):
 # plt.tight_layout()
 # plt.show()
 
+# function 1.1: regularly firing neuron
+def reg_neuron(n_timesteps, firing_rate, number=1):
+
+    # get the numner of seconds: n_timesteps is in ms
+    seconds = n_timesteps/1000
+
+    # get the number of spikes: firing rate is in Hz
+    n_spikes = seconds*firing_rate
+
+    # generate indices where spiking occurs
+    firing_indices = np.linspace(0,n_timesteps-1,n_spikes)
+
+    # make an array to reflect this activity
+    firing_neuron = np.zeros(n_timesteps)
+    firing_neuron[firing_indices]=1
+
+    print('number of spikes {}.'.format(np.size(np.nonzero(firing_neuron))))
+
+
+    return firing_neuron
+
+
+
+
+
 
 
 # Function 2: takes in an array of neuron activity and gives corresponding [NM]
@@ -127,22 +151,23 @@ def simulate_nm_conc(neuron_activity,nm_conc0, k_b,k_r,gamma):
 # plt.show()
 
 # Function 2.1: produces zoomed in plot
-def plot_nm_conc(nm,start,stop,colour='b', plotlabel = ''):
+def plot_nm_conc(nm,nmb, nmr, colour='b'):
 
     # define the timesteps to plot the [NM]
-    timesteps = stop - start + 1
-    t = np.linspace(start,stop,timesteps)
-
-    # get that section of the [NM] array
-    nm_section = nm[start:stop+1]
+    t = np.linspace(0,nm.size-1,nm.size)
 
     # plot the [NM] 
-    plt.plot(t,nm_section, color=colour, label=plotlabel)
+    plt.plot(t,nm, color=colour, label='[NM]')
+    plt.plot(t,nmb, color = 'g', label='[NM B]')
+    plt.plot(t,nmr, color = 'r', label='[NM R]')
     plt.xlabel('time (ms)')
     plt.ylabel('NM concentration')
-    plt.title('NM {} concentration from {} to {} ms'.format(plotlabel, start,stop))
+    plt.title('NM dynamics from firing neuron')
+    plt.legend()
     plt.show()
 
+
+
 
 
 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):
@@ -251,8 +276,8 @@ def exp_decay(t,a,b):
 
 
 
-
-def plot_f_signal(progression, progression_sub):
+# plot the progression from f -- df/f
+def plot_f_signal(progression, progression_sub, nm_conc):
 
 
     # create timesteps array for the plot
@@ -290,6 +315,39 @@ def plot_f_signal(progression, progression_sub):
     plt.suptitle('Progression from f to df/f', size = 16)
     plt.tight_layout()
 
+    plt.figure(2)
+    plt.subplot(2,2,1)
+    plt.plot(t,progression_sub[0], label='f')
+    plt.xlabel('time (ms)')
+    plt.ylabel('f')
+    plt.title('f  vs time')
+    plt.legend()
+
+    plt.subplot(2,2,2)
+    plt.plot(t,progression_sub[2], label = 'df')
+    plt.xlabel('time(ms)')
+    plt.ylabel(' df')
+    plt.title('df vs time (f0:median)')
+    plt.legend()
+
+    plt.subplot(2,2,3)
+    plt.plot(t,progression_sub[3], label = 'df/f')
+    plt.xlabel('time(ms)')
+    plt.ylabel(' df/f')
+    plt.title('df/f vs time (f0:median)')
+    plt.legend()
+
+    plt.subplot(2,2,4)
+    plt.plot(t,progression_sub[1], label='f - exponential (+ constant)')
+    plt.xlabel('time (ms)')
+    plt.ylabel('f_sub')
+    plt.title('f_sub vs time')
+    plt.legend()
+
+    plt.suptitle('Progression from f to df/f (using subtracted f)', size = 16)
+    plt.tight_layout()
+    plt.show()
+
 
 
 # # output 3