#=========================================================================== # Import, Set up MPI Variables, Load Necessary Files #=========================================================================== from mpi4py import MPI import time tic_0 = time.perf_counter() #script runtime calculation value import os from os.path import join import sys import numpy as np np.seterr(divide='ignore', invalid='ignore') from scipy import stats as st import neuron from neuron import h, gui import LFPy from LFPy import NetworkCell, Network, Synapse, RecExtElectrode, StimIntElectrode from matplotlib.collections import PolyCollection import matplotlib.pyplot as plt import pandas as pd #MPI variables: COMM = MPI.COMM_WORLD SIZE = COMM.Get_size() RANK = COMM.Get_rank() GLOBALSEED = int(sys.argv[1]) # Create new RandomState for each RANK SEED = GLOBALSEED*10000 np.random.seed(SEED + RANK) local_state = np.random.RandomState(SEED + RANK) halfnorm_rv = st.halfnorm halfnorm_rv.random_state = local_state uniform_rv = st.uniform uniform_rv.random_state = local_state #from net_params import * #Mechanisms and files print('Mechanisms found: ', os.path.isfile('mod/x86_64/special')) if RANK==0 else None neuron.h('forall delete_section()') neuron.load_mechanisms('mod/') h.load_file('net_functions.hoc') #=========================================================================== # Simulation, Analysis, and Plotting Controls #=========================================================================== TESTING = False # i.e.g generate 1 cell/pop, with 0.1 s runtime no_connectivity = False stimulate = False # Add a stimulus MDD = False #decrease PN GtonicApic and MN2PN weight by 40% DRUG = False rec_LFP = False #record LFP from center of layer rec_DIPOLES = False #record population - wide dipoles run_circuit_functions = True #=========================================================================== # Params #=========================================================================== dt = 0.025 #for both cell and network tstart = 0. tstop = 4500. celsius = 34. v_init = -80. #for both cell and network ###################### Load and organize Excel file ################################ circuit_params = {} #Import Excel file circuit_params = pd.read_excel('Circuit_param.xls', sheet_name = None, index_col = 0) #Get cell names and import biophys cell_names = [i for i in circuit_params['conn_probs'].axes[0]] for name in cell_names: h.load_file('models/biophys_'+name+'.hoc') circuit_params["syn_params"] = {'none':{'tau_r_AMPA': 0,'tau_d_AMPA': 0,'tau_r_NMDA': 0, 'tau_d_NMDA': 0, 'e': 0,'Dep': 0,'Fac': 0,'Use': 0,'u0':0,'gmax': 0}} circuit_params["multi_syns"] = {'none':{'loc':0,'scale':0}} # organizing dictionary for LFPY input for pre in cell_names: for post in cell_names: if "PYR" in pre: circuit_params["syn_params"][pre+post] = {'tau_r_AMPA': 0.3, 'tau_d_AMPA': 3, 'tau_r_NMDA': 2, 'tau_d_NMDA': 65, 'e': 0, 'u0':0, 'Dep': circuit_params["Depression"].at[pre, post], 'Fac': circuit_params["Facilitation"].at[pre, post], 'Use': circuit_params["Use"].at[pre, post], 'gmax': circuit_params["syn_cond"].at[pre, post]} else: circuit_params["syn_params"][pre+post] = {'tau_r': 1, 'tau_d': 10, 'e': -80, 'u0':0, 'Dep': circuit_params["Depression"].at[pre, post], 'Fac': circuit_params["Facilitation"].at[pre, post], 'Use': circuit_params["Use"].at[pre, post], 'gmax': circuit_params["syn_cond"].at[pre, post]} circuit_params["multi_syns"][pre+post] = {'loc':int(circuit_params["n_cont"].at[pre, post]),'scale':0} stimuli = [] for stimulus in circuit_params['STIM_PARAM'].axes[0]: stimuli.append({}) for param_name in circuit_params['STIM_PARAM'].axes[1]: stimuli[-1][param_name] = circuit_params['STIM_PARAM'].at[stimulus, param_name] new_param = circuit_params["syn_params"][stimuli[-1]['syn_params']].copy() new_param['gmax'] = stimuli[-1]['gmax'] stimuli[-1]['new_param'] = new_param COMM.Barrier() print('Importing, setting up MPI variables and loading necessary files took ', str((time.perf_counter() - tic_0)/60)[:5], 'minutes') if RANK==0 else None ################################################################################# if TESTING: OUTPUTPATH = 'Circuit_output_testing' for name in cell_names: circuit_params['SING_CELL_PARAM'].at['cell_num',name] = 1 tstop = 100 print('Running test...') if RANK ==0 else None else: OUTPUTPATH = 'Circuit_output' print('Running full simulation...') if RANK==0 else None COMM.Barrier() ################################################################################## networkParams = { 'dt' : dt, 'tstart': tstart, 'tstop' : tstop, 'v_init' : v_init, 'celsius' : celsius, 'OUTPUTPATH' : OUTPUTPATH, 'verbose': False} # L2/3 L4 L5 PYRmaxApics = [550 ,1550 ,1900] uppers = [-250 ,-1200 ,-1600] lowers = [-1200 ,-1580 ,-2300] depths = [] rangedepths = [] minSynLocs = [] syn_pos = [] pop_args = {} for i in range (3): depths.append((lowers[i]-uppers[i])/2-PYRmaxApics[i]) rangedepths.append(abs(lowers[i]-uppers[i])/2) minSynLocs.append((lowers[i]-uppers[i])/2*3-PYRmaxApics[i]) syn_pos.append({'section' : ['apic', 'dend'], 'fun' : [uniform_rv, halfnorm_rv], 'funargs' : [{'loc':minSynLocs[i], 'scale':abs(minSynLocs[i])},{'loc':minSynLocs[i], 'scale':abs(minSynLocs[i])}], 'funweights' : [1, 1.]}) syn_pos.append({'section' : ['apic'], 'fun' : [uniform_rv], 'funargs' : [{'loc':minSynLocs[i], 'scale':abs(minSynLocs[i])}], 'funweights' : [1.]}) syn_pos.append({'section' : ['dend'], 'fun' : [uniform_rv], 'funargs' : [{'loc':minSynLocs[i], 'scale':abs(minSynLocs[i])}], 'funweights' : [1.]}) syn_pos.append({'section' : ['dend'], 'fun' : [halfnorm_rv], 'funargs' : [{'loc':minSynLocs[i], 'scale':abs(minSynLocs[i])}], 'funweights' : [1.]}) names = ['HL2','HL4','HL5'] pop_args[names[i]]={'radius':250, 'loc':depths[i], 'scale':rangedepths[i]*4, 'cap':rangedepths[i]} # class RecExtElectrode parameters: L23_size = abs(uppers[1] - lowers[1]) e1 = 5 #-725 LFPelectrodeParameters = dict( x=np.zeros(1), y=np.zeros(1), z=[e1], N=np.array([[0., 1., 0.] for _ in range(1)]), r=5., n=50, sigma=0.3, method="soma_as_point") #method Network.simulate() parameters simargs = {'rec_imem': False, 'rec_vmem': False, 'rec_ipas': False, 'rec_icap': False, 'rec_isyn': False, 'rec_vmemsyn': False, 'rec_istim': False, 'rec_current_dipole_moment':rec_DIPOLES, 'rec_pop_contributions': False, 'rec_variables': [], 'to_memory': True, 'to_file': False, 'file_name':'OUTPUT.h5', 'dotprodcoeffs': None} #=========================================================================== # Functions #=========================================================================== def generateSubPop(popsize,mname,popargs,Gou,Gtonic,GtonicApic): print('Initiating ' + mname + ' population...') if RANK==0 else None morphpath = 'morphologies/' + mname + '.swc' templatepath = 'models/NeuronTemplate.hoc' templatename = 'NeuronTemplate' pt3d = True cellParams = { 'morphology': morphpath, 'templatefile': templatepath, 'templatename': templatename, 'templateargs': morphpath, 'v_init': v_init, #initial membrane potential, d=-65 'passive': False,#initialize passive mechs, d=T, should be overwritten by biophys 'dt': dt, 'tstart': 0., 'tstop': tstop,#defaults to 100 'nsegs_method': None, 'pt3d': pt3d,#use pt3d-info of the cell geometries switch, d=F 'delete_sections': False, 'verbose': False}#verbose output switch, for some reason doens't work, figure out why} rotation = {'x':circuit_params['SING_CELL_PARAM'].at['rotate_x', mname],'y':circuit_params['SING_CELL_PARAM'].at['rotate_y', mname]} popParams = { 'CWD': None, 'CELLPATH': None, 'Cell' : LFPy.NetworkCell, #play around with this, maybe put popargs into here 'POP_SIZE': int(popsize), 'name': mname, 'cell_args' : {**cellParams}, 'pop_args' : popargs, 'rotation_args' : rotation} network.create_population(**popParams) # Add biophys, OU processes, & tonic inhibition to cells for cellind in range(0,len(network.populations[mname].cells)): rseed = int(local_state.uniform()*SEED) biophys = 'h.biophys_' + mname + '(network.populations[\'' + mname + '\'].cells[' + str(cellind) + '].template)' exec(biophys) h.createArtificialSyn(rseed,network.populations[mname].cells[cellind].template,Gou) h.addTonicInhibition(network.populations[mname].cells[cellind].template,Gtonic,GtonicApic) def addStimulus(): cell_nums=[circuit_params['SING_CELL_PARAM'].at['cell_num',name] for name in cell_names] for stim in stimuli: stim_index = sum(cell_nums[:cell_names.index(stim['cell_name'])]) + stim['num_cells'] + stim['start_index'] for gid, cell in zip(network.populations[stim['cell_name']].gids, network.populations[stim['cell_name']].cells): if gid < stim_index and gid >= sum(cell_nums[:cell_names.index(stim['cell_name'])]) + stim['start_index']: idx = cell.get_rand_idx_area_norm(section=stim['loc'], nidx=stim['loc_num']) for i in idx: time_d=0 syn = Synapse(cell=cell, idx=i, syntype=stim['stim_type'], weight=1,**stim['new_param']) while time_d <= 0: time_d = np.random.uniform(low = stim['delay'], high = stim['delay']+stim['delay_range']) syn.set_spike_times_w_netstim(noise=0, start=(stim['start_time']+time_d), number=stim['num_stim'], interval=stim['interval'], seed=GLOBALSEED) #=========================================================================== # Sim #=========================================================================== network = Network(**networkParams) if MDD: # synaptic reduction for pre in cell_names: for post in cell_names: if 'SST' in pre: circuit_params["syn_params"][pre+post]["gmax"] = circuit_params["syn_cond"].at[pre, post]*0.6 # Synaptic reduction # tonic reduction for post in cell_names: if 'PYR' in post: circuit_params['SING_CELL_PARAM'].at['apic_tonic',post] = circuit_params['SING_CELL_PARAM'].at['apic_tonic',post]*0.6 circuit_params['SING_CELL_PARAM'].at["drug_apic_tonic",post] = circuit_params['SING_CELL_PARAM'].at["drug_apic_tonic",post]*0.6 else: sst = 0 total = 0 for pre in cell_names: if 'SST' in pre: sst += circuit_params["syn_cond"].at[pre, post]*circuit_params["n_cont"].at[pre,post]*circuit_params["conn_probs"].at[pre, post] total += circuit_params["syn_cond"].at[pre, post]*circuit_params["n_cont"].at[pre,post]*circuit_params["conn_probs"].at[pre, post] elif 'PV' in pre or 'VIP' in pre: total += circuit_params["syn_cond"].at[pre, post]*circuit_params["n_cont"].at[pre, post]*circuit_params["conn_probs"].at[pre, post] circuit_params['SING_CELL_PARAM'].at['norm_tonic',post] -= circuit_params['SING_CELL_PARAM'].at['norm_tonic',post]*sst/total*0.4 print(post, '_tonic reduced by: ', sst/total*100*0.4, '%') if RANK == 0 else None # Generate Populations tic = time.perf_counter() for cell_name in cell_names: if DRUG: generateSubPop(circuit_params['SING_CELL_PARAM'].at['cell_num',cell_name], cell_name,pop_args[cell_name[:3]], circuit_params['SING_CELL_PARAM'].at['GOU',cell_name], circuit_params['SING_CELL_PARAM'].at['drug_tonic',cell_name], circuit_params['SING_CELL_PARAM'].at['drug_apic_tonic',cell_name]) else: generateSubPop(circuit_params['SING_CELL_PARAM'].at['cell_num',cell_name], cell_name,pop_args[cell_name[:3]], circuit_params['SING_CELL_PARAM'].at['GOU',cell_name], circuit_params['SING_CELL_PARAM'].at['norm_tonic',cell_name], circuit_params['SING_CELL_PARAM'].at['apic_tonic',cell_name]) COMM.Barrier() print('Instantiating all populations took ', str((time.perf_counter() - tic_0)/60)[:5], 'minutes') if RANK==0 else None tic = time.perf_counter() # Synaptic Connection Parameters E_syn = neuron.h.ProbAMPANMDA I_syn = neuron.h.ProbUDFsyn for i, pre in enumerate(network.population_names): for j, post in enumerate(network.population_names): connectivity = network.get_connectivity_rand( pre=pre, post=post, connprob=0 if no_connectivity else circuit_params["conn_probs"].at[pre, post]) (conncount, syncount) = network.connect( pre=pre, post=post, connectivity=connectivity, syntype=E_syn if "PYR" in pre else I_syn, synparams=circuit_params["syn_params"][pre+post], weightfun=local_state.normal, weightargs={'loc':1, 'scale':0}, minweight=1, delayfun=local_state.normal, delayargs={'loc':0.5, 'scale':0}, mindelay=0.5, multapsefun=local_state.normal, multapseargs=circuit_params["multi_syns"][pre+post], syn_pos_args=syn_pos[circuit_params["Syn_pos"].at[pre,post]]) print('Connecting populations took ', str((time.perf_counter() - tic_0)/60)[:5], 'minutes') if RANK==0 else None # Setup Extracellular Recording Device COMM.Barrier() if stimulate: addStimulus() COMM.Barrier() # Run Simulation tic = time.perf_counter() if rec_LFP and not rec_DIPOLES: LFPelectrode = RecExtElectrode(**LFPelectrodeParameters) print('Simulating, recording SPIKES and LFP ... ') if RANK==0 else None SPIKES, OUTPUT, _ = network.simulate(electrode=LFPelectrode,**simargs) elif rec_LFP and rec_DIPOLES: print('Simulating, recording SPIKES, LFP, and DIPOLEMOMENTS ... ') if RANK==0 else None LFPelectrode = RecExtElectrode(**LFPelectrodeParameters) SPIKES, OUTPUT, DIPOLEMOMENT = network.simulate(electrode=LFPelectrode,**simargs) elif not rec_LFP and rec_DIPOLES: print('Simulating, recording SPIKES and DIPOLEMOMENTS ... ') if RANK==0 else None SPIKES, _, DIPOLEMOMENT = network.simulate(**simargs) elif not rec_LFP and not rec_DIPOLES: print('Simulating, recording SPIKES ... ') if RANK==0 else None SPIKES = network.simulate(**simargs) print('Simulation took ', str((time.perf_counter() - tic_0)/60)[:5], 'minutes') if RANK==0 else None COMM.Barrier() if RANK==0: tic = time.perf_counter() print('Saving simulation output ...') np.save(os.path.join(OUTPUTPATH,'SPIKES_Seed'+str(GLOBALSEED)+'.npy'), SPIKES) np.save(os.path.join(OUTPUTPATH,'OUTPUT_Seed'+str(GLOBALSEED)+'.npy'), OUTPUT) if rec_LFP else None np.save(os.path.join(OUTPUTPATH,'DIPOLEMOMENT_Seed'+str(GLOBALSEED)+'.npy'), DIPOLEMOMENT) if rec_DIPOLES else None print('Saving simulation took', str((time.perf_counter() - tic_0)/60)[:5], 'minutes') #=========================================================================== # Plotting #=========================================================================== if run_circuit_functions: tstart_plot = 2000 tstop_plot = tstop print('Creating/saving plots...') if RANK==0 else None exec(open("circuit_functions.py").read()) #=============== #Final printouts #=============== if TESTING: print('Test complete, switch TESTING to False for full simulation') if RANK==0 else None elif not TESTING: print('Simulation complete') if RANK==0 else None print('Script completed in ', str((time.perf_counter() - tic_0)/60)[:5], 'minutes') if RANK==0 else None