-*.png
-*.pkl
*.pyc
-params.json
out/
+venv
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.
+Terry Stewart, Xuan Choo, Chris Eliasmith. Some of those original files are
+preserved in old/.
<https://github.com/psipeter/drugs_and_working_memory>
+
+Usage
+=====
+
+Create and activate a new virtual environment:
+
+ python3 -m venv venv
+ . ./venv/bin/activate
+
+Install necesary dependencies:
+
+ pip install -r requirements.txt
+
+Modify conf.py as needed:
+
+ vi conf.py
+
+Run simulation (use the -h flag to see more options):
+
+ python3 model.py
+
+Results are stored in ./out/.
--- /dev/null
+import numpy as np
+import matplotlib.pyplot as plt
+from os import mkdir
+
+# 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):
+ try:
+ mkdir("./out")
+ except FileExistsError:
+ pass
+
+ out = f"./out/{self.__class__.__name__}"
+ 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(f"{out}-norep-activity.png", dpi=1000)
+
+ #######################################################################
+
+ plt.plot(self.x, self.gain)
+
+ plt.xlabel("Norepinephrine concentration (nM)")
+ plt.ylabel("Gain")
+ plt.title(f"Concentration vs Gain in {self.pretty}")
+
+ plt.draw()
+ plt.savefig(f"{out}-concentration-gain.png", dpi=1000)
+
+ #######################################################################
+
+ 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(f"{out}-concentration-bias.png", dpi=1000)
+
+class Alpha1(Alpha):
+ def __init__(self):
+ self.ki = 330
+ self.offset = 5.895
+ self.pretty = "α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 = "α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()
+++ /dev/null
-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()<P['misperceive']: perceived[n]=0
- return trials, perceived, cues
-
-
-'''drug approximations'''
-import nengo
-class MySolver(nengo.solvers.Solver):
- #When the simulator builds the network, it looks for a solver to calculate the decoders
- #instead of the normal least-squares solver, we define our own, so that we can return
- #the old decoders
- def __init__(self,weights): #feed in old decoders
- self.weights=False #they are not weights but decoders
- self.my_weights=weights
- def __call__(self,A,Y,rng=None,E=None): #the function that gets called by the builder
- return self.my_weights.T, dict()
-
-def reset_gain_bias(P,model,sim,wm,wm_recurrent,wm_to_decision,drug):
- #set gains and biases as a constant multiple of the old values
- wm.gain = sim.data[wm].gain * P['drug_effect_biophysical'][drug][0]
- wm.bias = sim.data[wm].bias * P['drug_effect_biophysical'][drug][1]
- #set the solver of each of the connections coming out of wm using the custom MySolver class
- #with input equal to the old decoders. We use the old decoders because we don't want the builder
- #to optimize the decoders to the new gainbias/bias values, otherwise it would "adapt" to the drug
- wm_recurrent.solver = MySolver(sim.model.params[wm_recurrent].weights)
- wm_to_decision.solver=MySolver(sim.model.params[wm_to_decision].weights)
- #rebuild the network to affect the gain/bias change
- sim=nengo.Simulator(model,dt=P['dt'])
- return sim
-
-'''dataframe initialization'''
-def primary_dataframe(P,sim,drug,trial,probe_wm,probe_output):
- import numpy as np
- import pandas as pd
- columns=('time','drug','wm','output','correct','trial')
- df_primary = pd.DataFrame(columns=columns, index=np.arange(0,len(P['timesteps'])))
- i=0
- for t in P['timesteps']:
- wm_val=np.abs(sim.data[probe_wm][t][0])
- output_val=sim.data[probe_output][t][0]
- correct=get_correct(P['cues'][trial],output_val)
- rt=t*P['dt_sample']
- df_primary.loc[i]=[rt,drug,wm_val,output_val,correct,trial]
- i+=1
- return df_primary
-
-def firing_dataframe(P,sim,drug,trial,sim_wm,probe_spikes):
- import numpy as np
- import pandas as pd
- columns=('time','drug','neuron-trial','tuning','firing_rate')
- df_firing = pd.DataFrame(columns=columns, index=np.arange(0,len(P['timesteps'])*\
- int(P['neurons_wm']*P['frac'])))
- t_width = 0.2
- t_h = np.arange(t_width / P['dt']) * P['dt'] - t_width / 2.0
- h = np.exp(-t_h ** 2 / (2 * P['sigma_smoothing'] ** 2))
- h = h / np.linalg.norm(h, 1)
- j=0
- for nrn in range(int(P['neurons_wm']*P['frac'])):
- enc = sim_wm.encoders[nrn]
- tuning = get_tuning(P,trial,enc)
- spikes = sim.data[probe_spikes][:,nrn]
- firing_rate = np.convolve(spikes,h,mode='same')
- for t in P['timesteps']:
- rt=t*P['dt_sample']
- df_firing.loc[j]=[rt,drug,nrn+trial*P['neurons_wm'],tuning,firing_rate[t]]
- j+=1
- # print 'appending dataframe for neuron %s' %f
- return df_firing
-
-def get_correct(cue,output_val):
- if (cue > 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
+++ /dev/null
-'''
-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
+++ /dev/null
-import numpy as np
-import matplotlib.pyplot as plt
-from os import mkdir
-
-# 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):
- try:
- mkdir("./out")
- except FileExistsError:
- pass
-
- out = f"./out/{self.__class__.__name__}"
- 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(f"{out}-norep-activity.png", dpi=1000)
-
- #######################################################################
-
- plt.plot(self.x, self.gain)
-
- plt.xlabel("Norepinephrine concentration (nM)")
- plt.ylabel("Gain")
- plt.title(f"Concentration vs Gain in {self.pretty}")
-
- plt.draw()
- plt.savefig(f"{out}-concentration-gain.png", dpi=1000)
-
- #######################################################################
-
- 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(f"{out}-concentration-bias.png", dpi=1000)
-
-class Alpha1(Alpha):
- def __init__(self):
- self.ki = 330
- self.offset = 5.895
- self.pretty = "α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 = "α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()
+++ /dev/null
-{
-'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
--- /dev/null
+*
+!.py
+!parameters.txt
--- /dev/null
+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()<P['misperceive']: perceived[n]=0
+ return trials, perceived, cues
+
+
+'''drug approximations'''
+import nengo
+class MySolver(nengo.solvers.Solver):
+ #When the simulator builds the network, it looks for a solver to calculate the decoders
+ #instead of the normal least-squares solver, we define our own, so that we can return
+ #the old decoders
+ def __init__(self,weights): #feed in old decoders
+ self.weights=False #they are not weights but decoders
+ self.my_weights=weights
+ def __call__(self,A,Y,rng=None,E=None): #the function that gets called by the builder
+ return self.my_weights.T, dict()
+
+def reset_gain_bias(P,model,sim,wm,wm_recurrent,wm_to_decision,drug):
+ #set gains and biases as a constant multiple of the old values
+ wm.gain = sim.data[wm].gain * P['drug_effect_biophysical'][drug][0]
+ wm.bias = sim.data[wm].bias * P['drug_effect_biophysical'][drug][1]
+ #set the solver of each of the connections coming out of wm using the custom MySolver class
+ #with input equal to the old decoders. We use the old decoders because we don't want the builder
+ #to optimize the decoders to the new gainbias/bias values, otherwise it would "adapt" to the drug
+ wm_recurrent.solver = MySolver(sim.model.params[wm_recurrent].weights)
+ wm_to_decision.solver=MySolver(sim.model.params[wm_to_decision].weights)
+ #rebuild the network to affect the gain/bias change
+ sim=nengo.Simulator(model,dt=P['dt'])
+ return sim
+
+'''dataframe initialization'''
+def primary_dataframe(P,sim,drug,trial,probe_wm,probe_output):
+ import numpy as np
+ import pandas as pd
+ columns=('time','drug','wm','output','correct','trial')
+ df_primary = pd.DataFrame(columns=columns, index=np.arange(0,len(P['timesteps'])))
+ i=0
+ for t in P['timesteps']:
+ wm_val=np.abs(sim.data[probe_wm][t][0])
+ output_val=sim.data[probe_output][t][0]
+ correct=get_correct(P['cues'][trial],output_val)
+ rt=t*P['dt_sample']
+ df_primary.loc[i]=[rt,drug,wm_val,output_val,correct,trial]
+ i+=1
+ return df_primary
+
+def firing_dataframe(P,sim,drug,trial,sim_wm,probe_spikes):
+ import numpy as np
+ import pandas as pd
+ columns=('time','drug','neuron-trial','tuning','firing_rate')
+ df_firing = pd.DataFrame(columns=columns, index=np.arange(0,len(P['timesteps'])*\
+ int(P['neurons_wm']*P['frac'])))
+ t_width = 0.2
+ t_h = np.arange(t_width / P['dt']) * P['dt'] - t_width / 2.0
+ h = np.exp(-t_h ** 2 / (2 * P['sigma_smoothing'] ** 2))
+ h = h / np.linalg.norm(h, 1)
+ j=0
+ for nrn in range(int(P['neurons_wm']*P['frac'])):
+ enc = sim_wm.encoders[nrn]
+ tuning = get_tuning(P,trial,enc)
+ spikes = sim.data[probe_spikes][:,nrn]
+ firing_rate = np.convolve(spikes,h,mode='same')
+ for t in P['timesteps']:
+ rt=t*P['dt_sample']
+ df_firing.loc[j]=[rt,drug,nrn+trial*P['neurons_wm'],tuning,firing_rate[t]]
+ j+=1
+ # print 'appending dataframe for neuron %s' %f
+ return df_firing
+
+def get_correct(cue,output_val):
+ if (cue > 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
--- /dev/null
+'''
+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
--- /dev/null
+{
+'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
--- /dev/null
+autopep8==1.5.7
+cycler==0.10.0
+kiwisolver==1.1.0
+matplotlib==3.4.2
+nengo==3.1.0
+numpy==1.21.1
+pandas==1.3.1
+Pillow==8.3.1
+pycodestyle==2.7.0
+pyparsing==2.4.7
+python-dateutil==2.8.2
+pytz==2021.1
+scipy==1.5.3
+six==1.16.0
+toml==0.10.2
projects / norepinephrine_wm.git / commitdiff