From 3a97d6c3ff2de0dac3add612961b2f6abb273531 Mon Sep 17 00:00:00 2001 From: Armaan Bhojwani Date: Fri, 6 Aug 2021 18:42:03 -0400 Subject: [PATCH] Initial commit --- .gitignore | 5 + LICENSE | 1 + README | 12 ++ norepinephrine_wm/helper.py | 131 +++++++++++++++++++ norepinephrine_wm/model.py | 208 +++++++++++++++++++++++++++++++ norepinephrine_wm/newmodel.py | 111 +++++++++++++++++ norepinephrine_wm/parameters.txt | 32 +++++ 7 files changed, 500 insertions(+) create mode 100644 .gitignore create mode 100644 LICENSE create mode 100644 README create mode 100644 norepinephrine_wm/helper.py create mode 100644 norepinephrine_wm/model.py create mode 100644 norepinephrine_wm/newmodel.py create mode 100644 norepinephrine_wm/parameters.txt diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..0ef23ea --- /dev/null +++ b/.gitignore @@ -0,0 +1,5 @@ +*.png +*.pkl +*.pyc +params.json +out/ diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..55bae75 --- /dev/null +++ b/LICENSE @@ -0,0 +1 @@ +This codebase was originally unlicensed. diff --git a/README b/README new file mode 100644 index 0000000..9b1367a --- /dev/null +++ b/README @@ -0,0 +1,12 @@ +Modeling the Effects of Norepinephrine Concentration on Working Memory in the +Prefrontal Cortex with a1 and a2 Receptors. + +Python model doing the above. Part of Boston University RISE 2021. + +Julia Tomaszewska, Armaan Bhojwani, Thalia Inui, Nitsueh Kebere, Dr. +Marianne Bezaire. + +The work in this repository was in-part derived from "Effects of Guanfacine and +Phenylephrine on a Spiking Neuron Model of Working Memory" by Peter Duggins, +Terry Stewart, Xuan Choo, Chris Eliasmith. + diff --git a/norepinephrine_wm/helper.py b/norepinephrine_wm/helper.py new file mode 100644 index 0000000..078b2f1 --- /dev/null +++ b/norepinephrine_wm/helper.py @@ -0,0 +1,131 @@ +def make_addon(N): + import string + import random + addon=str(''.join(random.choice(string.ascii_uppercase+string.digits) for _ in range(N))) + return addon + +def ch_dir(): + #change directory for data and plot outputs + import os + import sys + root=os.getcwd() + addon=make_addon(9) + datadir='' + if sys.platform == "linux" or sys.platform == "linux2" or sys.platform == "darwin": + datadir=root+'/data/'+addon+'/' #linux or mac + elif sys.platform == "win32": + datadir=root+'\\data\\'+addon+'\\' #windows + os.makedirs(datadir) + os.chdir(datadir) + return datadir + +def empirical_dataframe(): + import numpy as np + import pandas as pd + columns=('time','drug','accuracy','trial') + emp_times = [3.0,5.0,7.0,9.0] + emp_dataframe = pd.DataFrame(columns=columns,index=np.arange(0, 12)) + pre_PHE=[0.972, 0.947, 0.913, 0.798] + pre_GFC=[0.970, 0.942, 0.882, 0.766] + post_GFC=[0.966, 0.928, 0.906, 0.838] + post_PHE=[0.972, 0.938, 0.847, 0.666] + q=0 + for t in range(len(emp_times)): + emp_dataframe.loc[q]=[emp_times[t],'control (empirical)',np.average([pre_GFC[t],pre_PHE[t]]),0] + emp_dataframe.loc[q+1]=[emp_times[t],'PHE (empirical)',post_PHE[t],0] + emp_dataframe.loc[q+2]=[emp_times[t],'GFC (empirical)',post_GFC[t],0] + q+=3 + return emp_dataframe + +def make_cues(P): + import numpy as np + trials=np.arange(P['n_trials']) + perceived=np.ones(P['n_trials']) #list of correctly perceived (not necessarily remembered) cues + rng=np.random.RandomState(seed=P['seed']) + cues=2*rng.randint(2,size=P['n_trials'])-1 #whether the cues is on the left or right + for n in range(len(perceived)): + if rng.rand() 0.0 and output_val > 0.0) or (cue < 0.0 and output_val < 0.0): correct=1 + else: correct=0 + return correct + +def get_tuning(P,trial,enc): + cue=P['cues'][trial] + enc_min_cutoff=P['enc_min_cutoff'] + enc_max_cutoff=P['enc_max_cutoff'] + if (cue > 0.0 and 0.0 < enc[0] < enc_min_cutoff) or \ + (cue < 0.0 and 0.0 > enc[0] > -1.0*enc_min_cutoff): tuning='superweak' + if (cue > 0.0 and enc_min_cutoff < enc[0] < enc_max_cutoff) or \ + (cue < 0.0 and -1.0*enc_max_cutoff < enc[0] < -1.0*enc_min_cutoff): tuning='weak' + elif (cue > 0.0 and enc[0] > enc_max_cutoff) or \ + (cue < 0.0 and enc[0] < -1.0*enc_max_cutoff): tuning='strong' + else: tuning='nonpreferred' + return tuning \ No newline at end of file diff --git a/norepinephrine_wm/model.py b/norepinephrine_wm/model.py new file mode 100644 index 0000000..a81eb76 --- /dev/null +++ b/norepinephrine_wm/model.py @@ -0,0 +1,208 @@ +''' +Peter Duggins, Terry Stewart, Xuan Choo, Chris Eliasmith +Effects of Guanfacine and Phenylephrine on a Spiking Neuron Model of Working Memory +June-August 2016 +Main Model File +''' + +def run(params): + import nengo + from nengo.rc import rc + import numpy as np + import pandas as pd + from helper import reset_gain_bias, primary_dataframe, firing_dataframe + + decision_type=params[0] + drug_type=params[1] + drug = params[2] + trial = params[3] + seed = params[4] + P = params[5] + dt=P['dt'] + dt_sample=P['dt_sample'] + t_cue=P['t_cue'] + t_delay=P['t_delay'] + drug_effect_neural=P['drug_effect_neural'] + drug_effect_functional=P['drug_effect_functional'] + drug_effect_biophysical=P['drug_effect_biophysical'] + enc_min_cutoff=P['enc_min_cutoff'] + enc_max_cutoff=P['enc_max_cutoff'] + sigma_smoothing=P['sigma_smoothing'] + frac=P['frac'] + neurons_inputs=P['neurons_inputs'] + neurons_wm=P['neurons_wm'] + neurons_decide=P['neurons_decide'] + time_scale=P['time_scale'] + cue_scale=P['cue_scale'] + tau=P['tau'] + tau_wm=P['tau_wm'] + noise_wm=P['noise_wm'] + noise_decision=P['noise_decision'] + perceived=P['perceived'] + cues=P['cues'] + + if drug_type == 'biophysical': rc.set("decoder_cache", "enabled", "False") #don't try to remember old decoders + else: rc.set("decoder_cache", "enabled", "True") + + def cue_function(t): + if t < t_cue and perceived[trial]!=0: + return cue_scale * cues[trial] + else: return 0 + + def time_function(t): + if t > t_cue: + return time_scale + else: return 0 + + def noise_bias_function(t): + import numpy as np + if drug_type=='neural': + return np.random.normal(drug_effect_neural[drug],noise_wm) + else: + return np.random.normal(0.0,noise_wm) + + def noise_decision_function(t): + import numpy as np + if decision_type == 'default': + return np.random.normal(0.0,noise_decision) + elif decision_type == 'basal_ganglia': + return np.random.normal(0.0,noise_decision,size=2) + + def inputs_function(x): + return x * tau_wm + + def wm_recurrent_function(x): + if drug_type == 'functional': + return x * drug_effect_functional[drug] + else: + return x + + def decision_function(x): + output=0.0 + if decision_type=='default': + value=x[0]+x[1] + if value > 0.0: output = 1.0 + elif value < 0.0: output = -1.0 + elif decision_type=='basal_ganglia': + if x[0] > x[1]: output = 1.0 + elif x[0] < x[1]: output = -1.0 + return output + + def BG_rescale(x): #rescales -1 to 1 into 0.3 to 1, makes 2-dimensional + pos_x = 0.5 * (x + 1) + rescaled = 0.4 + 0.6 * pos_x, 0.4 + 0.6 * (1 - pos_x) + return rescaled + + '''model definition''' + with nengo.Network(seed=seed+trial) as model: + + #Ensembles + cue = nengo.Node(output=cue_function) + time = nengo.Node(output=time_function) + inputs = nengo.Ensemble(neurons_inputs,2) + noise_wm_node = nengo.Node(output=noise_bias_function) + noise_decision_node = nengo.Node(output=noise_decision_function) + wm = nengo.Ensemble(neurons_wm,2) + decision = nengo.Ensemble(neurons_decide,2) + output = nengo.Ensemble(neurons_decide,1) + + #Connections + nengo.Connection(cue,inputs[0],synapse=None) + nengo.Connection(time,inputs[1],synapse=None) + nengo.Connection(inputs,wm,synapse=tau_wm,function=inputs_function) + wm_recurrent=nengo.Connection(wm,wm,synapse=tau_wm,function=wm_recurrent_function) + nengo.Connection(noise_wm_node,wm.neurons,synapse=tau_wm,transform=np.ones((neurons_wm,1))*tau_wm) + wm_to_decision=nengo.Connection(wm[0],decision[0],synapse=tau) + nengo.Connection(noise_decision_node,decision[1],synapse=None) + nengo.Connection(decision,output,function=decision_function) + + #Probes + probe_wm=nengo.Probe(wm[0],synapse=0.01,sample_every=dt_sample) + probe_spikes=nengo.Probe(wm.neurons, 'spikes', sample_every=dt_sample) + probe_output=nengo.Probe(output,synapse=None,sample_every=dt_sample) + + + + + '''SIMULATION''' + print 'Running drug \"%s\", trial %s...' %(drug,trial+1) + with nengo.Simulator(model,dt=dt) as sim: + if drug_type == 'biophysical': sim=reset_gain_bias( + P,model,sim,wm,wm_recurrent,wm_to_decision,drug) + sim.run(t_cue+t_delay) + df_primary=primary_dataframe(P,sim,drug,trial,probe_wm,probe_output) + df_firing=firing_dataframe(P,sim,drug,trial,sim.data[wm],probe_spikes) + return [df_primary, df_firing] + + + +'''MAIN''' +def main(): + import matplotlib.pyplot as plt + import seaborn as sns + import pandas as pd + import numpy as np + from helper import make_cues, empirical_dataframe + from pathos.helpers import freeze_support #for Windows + # import ipdb + + '''Import Parameters from File''' + P=eval(open('parameters.txt').read()) #parameter dictionary + seed=P['seed'] #sets tuning curves equal to control before drug application + n_trials=P['n_trials'] + drug_type=str(P['drug_type']) + decision_type=str(P['decision_type']) + drugs=P['drugs'] + trials, perceived, cues = make_cues(P) + P['timesteps']=np.arange(0,int((P['t_cue']+P['t_delay'])/P['dt_sample'])) + P['cues']=cues + P['perceived']=perceived + + '''Multiprocessing''' + print "drug_type=%s, decision_type=%s, trials=%s..." %(drug_type,decision_type,n_trials) + freeze_support() + exp_params=[] + for drug in drugs: + for trial in trials: + exp_params.append([decision_type, drug_type, drug, trial, seed, P]) + df_list=[run(exp_params[0]),run(exp_params[-1])] + primary_dataframe = pd.concat([df_list[i][0] for i in range(len(df_list))], ignore_index=True) + firing_dataframe = pd.concat([df_list[i][1] for i in range(len(df_list))], ignore_index=True) + + '''Plot and Export''' + print 'Exporting Data...' + primary_dataframe.to_pickle('primary_data.pkl') + firing_dataframe.to_pickle('firing_data.pkl') + param_df=pd.DataFrame([P]) + param_df.reset_index().to_json('params.json',orient='records') + + print 'Plotting...' + emp_dataframe=empirical_dataframe() + sns.set(context='poster') + figure, (ax1, ax2) = plt.subplots(2, 1) + sns.tsplot(time="time",value="wm",data=primary_dataframe,unit="trial",condition='drug',ax=ax1,ci=95) + sns.tsplot(time="time",value="correct",data=primary_dataframe,unit="trial",condition='drug',ax=ax2,ci=95) + sns.tsplot(time="time",value="accuracy",data=emp_dataframe,unit='trial',condition='drug', + interpolate=False,ax=ax2) + sns.tsplot(time="time",value="accuracy",data=emp_dataframe, unit='trial',condition='drug', + interpolate=True,ax=ax2, legend=False) + ax1.set(xlabel='',ylabel='decoded $\hat{cue}$',xlim=(0,9.5),ylim=(0,1), + title="drug_type=%s, decision_type=%s, trials=%s" %(drug_type,decision_type,n_trials)) + ax2.set(xlabel='time (s)',xlim=(0,9.5),ylim=(0.5,1),ylabel='DRT accuracy') + figure.savefig('primary_plots.png') + + figure2, (ax3, ax4) = plt.subplots(1, 2) + if len(firing_dataframe.query("tuning=='weak'"))>0: + sns.tsplot(time="time",value="firing_rate",unit="neuron-trial",condition='drug',ax=ax3,ci=95, + data=firing_dataframe.query("tuning=='weak'").reset_index()) + if len(firing_dataframe.query("tuning=='nonpreferred'"))>0: + sns.tsplot(time="time",value="firing_rate",unit="neuron-trial",condition='drug',ax=ax4,ci=95, + data=firing_dataframe.query("tuning=='nonpreferred'").reset_index()) + ax3.set(xlabel='time (s)',xlim=(0.0,9.5),ylim=(0,250),ylabel='Normalized Firing Rate',title='Preferred Direction') + ax4.set(xlabel='time (s)',xlim=(0.0,9.5),ylim=(0,250),ylabel='',title='Nonpreferred Direction') + figure2.savefig('firing_plots.png') + + plt.show() + +if __name__=='__main__': + main() \ No newline at end of file diff --git a/norepinephrine_wm/newmodel.py b/norepinephrine_wm/newmodel.py new file mode 100644 index 0000000..c6afe03 --- /dev/null +++ b/norepinephrine_wm/newmodel.py @@ -0,0 +1,111 @@ +#!/usr/bin/env python2 +# -*- coding: utf-8 -*- + +import numpy as np +import matplotlib.pyplot as plt +import nengo + +# Refer to parameters.txt for documentation +dt = 0.001 +t_cue = 1.0 +cue_scale = 1.0 +perceived = 0 # ??? +time_scale = 0.4 +steps = 100 + +class Alpha(object): + def __init__(self): + self.x = np.logspace(0, 3, steps) + self.y = 1 / (1 + (999 * np.exp(-0.1233 * (self.x / self.offset)))) + + self.gain = [] + self.bias = [] + + def calcgb(self, gaind, biasd): + for i in range(steps): + y = self.y[i] + self.gain.append(1 + gaind * y) + self.bias.append(1 + biasd * y) + + def plot(self): + plt.plot(self.x, self.y) + + plt.xlabel("Norepinephrine concentration (nM)") + plt.ylabel("Activity (%)") + plt.title("Norepinepherine Concentration vs Neuron Activity in " + self.pretty) + + plt.vlines(self.ki, 0, 1, linestyles="dashed") + plt.text(1.1 * self.ki, 0.1, "Affinity") + + plt.hlines(0.5, 0, 1000, linestyles="dashed") + plt.text(1, 0.51, "50%") + + plt.xscale("log") + gc = plt.gca() + gc.set_yticklabels(['{:.0f}%'.format(x * 100) for x in gc.get_yticks()]) + + plt.draw() + plt.savefig(self.__class__.__name__ + "-norep-activity.png", dpi=1000) + plt.show() + + ####################################################################### + + plt.plot(self.x, self.gain) + + plt.xlabel("Norepinephrine concentration (nM)") + plt.ylabel("Gain") + plt.title("Concentration vs Gain in " + self.pretty) + + plt.draw() + plt.savefig(self.__class__.__name__ + "-concentration-gain.png", dpi=1000) + plt.show() + + ####################################################################### + + plt.plot(self.x, self.bias) + + plt.xlabel("Norepinephrine concentration (nM)") + plt.ylabel("Bias") + plt.title("Concentration vs Bias in " + self.pretty) + + plt.draw() + plt.savefig(self.__class__.__name__ + "-concentration-bias.png", dpi=1000) + plt.show() + +class Alpha1(Alpha): + def __init__(self): + self.ki = 330 + self.offset = 5.895 + self.pretty = u"α1 Receptor" + self.gaind = 0.1 + self.biasd = 0.1 + super().__init__() + + def calcgb(self): + super().calcgb(self.gaind, self.biasd) + +class Alpha2(Alpha): + def __init__(self): + self.ki = 56 + self.offset = 1 + self.pretty = u"α2 Receptor" + self.gaind = -0.04 + self.biasd = -0.02 + super().__init__() + + def calcgb(self): + super().calcgb(self.gaind, self.biasd) + +def main(): + plt.style.use("ggplot") + + a1 = Alpha1() + a1.calcgb() + a1.plot() + + a2 = Alpha2() + a2.calcgb() + a2.plot() + +if __name__ == "__main__": + main() diff --git a/norepinephrine_wm/parameters.txt b/norepinephrine_wm/parameters.txt new file mode 100644 index 0000000..d36549b --- /dev/null +++ b/norepinephrine_wm/parameters.txt @@ -0,0 +1,32 @@ +{ +'seed':3, #for the simulator build process +'n_trials':3, +'n_processes':9, +'dt':0.001, #timestep +'dt_sample':0.01, #timestep for data recording through probes +'t_cue':1.0, #duration of cue presentation +'t_delay':8.0, #duration of delay period between cue and decision +'misperceive':0.1, #chance of failing to perceive the cue, causing no info to go into WM + +'neurons_inputs':100, #neurons for the inputs ensemble +'neurons_wm':1000, #neurons for workimg memory ensemble +'neurons_decide':100, #neurons for decision or basal ganglia +'time_scale':0.4, #how fast does the 'time' dimension accumulate in WM neurons +'cue_scale':1.0, #how strong is the cueulus from the visual system +'tau':0.01, #synaptic time constant between ensembles +'tau_wm':0.1, #synapse on recurrent connection in wm +'noise_wm':0.005, #std of full-spectrum white noise inject into wm.neurons +'noise_decision':0.3, #std of added gaussian noise in the default decision procedure; + +'decision_type':'default', #decision procedure: 'default', 'basal_ganglia' +'drug_type':'biophysical', #drug simulation: 'neural','functional','biophysical' +'drugs':['control','PHE','GFC'], #names of the of drugs to simulate, see dicts below +'drug_effect_neural':{'control':0.0,'PHE':-0.2,'GFC':0.5}, #mean of cue onto wm.neurons +'drug_effect_biophysical':{'control':[1.0,1,0],'PHE':[0.99,1.02],'GFC':[1.05,0.95]}, #k_gain, k_bias +'drug_effect_functional':{'control':1.0,'PHE':0.985,'GFC':1.03}, #k_multiply + +'enc_min_cutoff':0.3, #minimum cutoff for "weak" encoders in preferred directions +'enc_max_cutoff':0.6, #maximum cutoff for "weak" encoders in preferred directions +'sigma_smoothing':0.005, #gaussian smoothing applied to spike data to calculate firing rate +'frac':0.01, #fraction of neurons in WM to add to firing rate dataframe and plot +} \ No newline at end of file -- 2.39.2