#This script is used to carry out neuron simulations described in 

#Knox AT, Thompson CH, Scott D, Abramova TV, Stieve B, Freeman A, George AL Jr. Genotype-function-phenotype correlations for SCN1A variants identified by clinical genetic testing. Ann Clin Transl Neurol. 2025 Mar;12(3):499-511. doi: 10.1002/acn3.52297. Epub 2025 Jan 21. PMID: 39838578; PMCID: PMC11920720.

#The original PV+ Interneuron model on which this model is based can be found at https://modeldb.science/264834 and is described in the publication

#Berecki G, Bryson A, Terhag J, Maljevic S, Gazina EV, Hill SL, Petrou S. SCN1A gain of function in early infantile encephalopathy. Ann Neurol. 2019 Apr;85(4):514-525. doi: 10.1002/ana.25438. Epub 2019 Mar 7. PMID: 30779207.

#Written by Andrew Knox

from netpyne import specs, sim
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import pandas as pd
import random
import sys
import os
import seaborn as sns
from scipy import stats

#todo: fix sloppy interchangable use of g and gid

if len(sys.argv) > 1:
    seed = sys.argv[1]
else:
    seed = 1 #default seed

eiratio = 3.0 
epsp_rate = 40.0 
sec_to_stim = 'soma_0'
sim_duration = 200.0
record_sec = 'soma_0'

putInTopDirectory = True 
stimulation_paradigm = 'epspipspsoma' #options iclampsoma, epspsoma, epspipspsoma,epspipspdendrites

netParams = specs.NetParams()
simConfig = specs.SimConfig()

#setup neurons
netParams.popParams['WT'] = {
    "cellType": "INT-WT",
    "numCells":1,
    "xnormRange":[0,1],
    "znormRange":[0,1],
    "cellModel":"INT3D"
}

netParams.popParams['I1356M'] = {
    "cellType": "INT-I1356M",
    "numCells": 1,
    "xnormRange":[1.5,2.5],
    "znormRange":[1.5,2.5],
    "cellModel":"INT3D"
}

netParams.popParams['L479P'] = {
    "cellType": "INT-L479P",
    "numCells": 1,
    "xnormRange":[4,5],
    "znormRange":[4,5],
    "cellModel":"INT3D"
}

netParams.popParams['KO'] = {
    "cellType": "INT-KO",
    "numCells": 1,
    "xnormRange":[7,8],
    "znormRange":[7,8],
    "cellModel":"INT3D"
}

netParams.importCellParams(
        label='INT-WT', #previously interneuron_hoc
        conds={'cellType': 'INT-WT','cellModel':'INT3D'},
        fileName='PV_interneuron.hoc',
        cellName="BC")

netParams.importCellParams(
        label='INT-I1356M', #'interneuron_I1356_hoc', 
        conds={'cellType': 'INT-I1356M','cellModel':'INT3D'},
        fileName='PV_interneuron.hoc',
        cellName="BC")

netParams.importCellParams(
        label='INT-L479P',   #'interneuron_L479P_hoc',
        conds={'cellType': 'INT-L479P','cellModel':'INT3D'},
        fileName='PV_interneuron.hoc',
        cellName="BC")

netParams.importCellParams(
        label='INT-KO', #'interneuron_KO_hoc',
        conds={'cellType': 'INT-KO','cellModel':'INT3D'},
        fileName='PV_interneuron.hoc',
        cellName="BC")

netParams.defaultThreshold = 0.0

netParams.synMechParams['AMPA'] = {
    "mod": "Exp2Syn",
    "tau1": 0.5,
    "tau2": 2.4,
    "e": 0
}
netParams.synMechParams['GABA'] = {
    "mod": "Exp2Syn",
    "tau1": 1.0,
    "tau2": 7.0,
    "e": -70
}


#setup network
simConfig.duration = sim_duration 
simConfig.dt = 0.1
simConfig.hParams = {
    "celsius": 6.3,
    "v_init": -65.0,
    "clamp_resist": 0.001
}
simConfig.recordCells = [
    0,
    1,
    2,
    3
]
simConfig.recordTraces = {
    "V_soma": {
        "sec": record_sec,
        "loc": 0.5,
        "var": "v"
    }
}

simConfig.recordStim = True

#create network
sim.create(netParams = netParams, simConfig = simConfig)
record_gids = sim.cfg.recordCells

sectionlist = sim.net.cells[0].secs.keys()

dendlist = [s for s in sectionlist if "dend" in s]
somalist = [s for s in sectionlist if "soma" in s]
axonlist = [s for s in sectionlist if "axon" in s]

epspid='0'
ipspid='0'
simid='0'

secs_to_stim_epsp = random.sample(dendlist,6)
secs_to_stim_ipsp = random.sample(dendlist,6)

all_stim_secs = secs_to_stim_epsp + secs_to_stim_ipsp

# code to determine whether excitatory or inhibitory synapses are closer to soma, or whether it's a mixed picture
def updateNodeType(neuron, node, child, node_status_dict, all_stim_secs):
    #print("node:",node)
    if node in all_stim_secs or node == 'soma_0':
        return node_status_dict 

    n_status = node_status_dict[node]
    c_status = node_status_dict[child]
    
    if n_status == 'M':
        return node_status_dict
    elif n_status == 'none':
        node_status_dict[node] = c_status
    elif (n_status == 'I' and c_status == 'E') or (n_status == 'E' and c_status == 'I') or c_status == 'M':
        node_status_dict[node] = 'M'

    parent = neuron.secs[node].topol['parentSec']
    #print("call", parent, node, node_status_dict)
    updateNodeType(neuron, parent, node, node_status_dict, all_stim_secs)
    return node_status_dict

node_status_dict = {}
for d in dendlist:
    node_status_dict[d] = 'none'

for s in secs_to_stim_epsp:
    node_status_dict[s] = 'E'

for s in secs_to_stim_ipsp:
    node_status_dict[s] = 'I'

for sec in all_stim_secs:
    parent = sim.net.cells[0].secs[sec].topol['parentSec']
    print(parent)
    print("nsd:",node_status_dict)
    node_status_dict = updateNodeType(sim.net.cells[0],parent,sec,node_status_dict,all_stim_secs)
    print("nsda:",node_status_dict)

top_sec_list = []

for sec in dendlist:
    if sim.net.cells[0].secs[sec].topol['parentSec'] == 'soma_0':
        top_sec_list.append(sec)

top_sec_exc_list = []

syn_exc_rating = 0
for sec in top_sec_list:
    if node_status_dict[sec] == 'E':
        top_sec_exc_list.append((sec,'E'))
        syn_exc_rating += 1
    elif node_status_dict[sec] == 'I':
        top_sec_exc_list.append((sec,'I'))
        syn_exc_rating -= 1
    elif node_status_dict[sec] == 'M':
        top_sec_exc_list.append((sec,'M'))
        syn_exc_rating += 0.5
    else:
        top_sec_exc_list.append((sec,'none'))

#print("top:",top_sec_exc_list)
#print("Rating:",syn_exc_rating)
syn_dist_rating = 0

for sec in secs_to_stim_epsp:
    parent = sim.net.cells[0].secs[sec].topol['parentSec']
    parent_list = []
    parent_list.append(parent)
    while parent != 'soma_0':
        parent = sim.net.cells[0].secs[parent].topol['parentSec']
        parent_list.append(parent)
    syn_dist_rating += len(parent_list)

for sec in secs_to_stim_ipsp:
    parent = sim.net.cells[0].secs[sec].topol['parentSec']
    parent_list = []
    parent_list.append(parent)
    while parent != 'soma_0':
        parent = sim.net.cells[0].secs[parent].topol['parentSec']
        parent_list.append(parent)
    syn_dist_rating -= len(parent_list)


#set up cell stimulation based on the chosen paradigm (iclamp vs epsp vs ipsp, dendrites vs soma
for g,gid in enumerate(record_gids):
    cellPop = sim.net.cells[gid].tags['pop']
    print("pop:",cellPop)
    cellType = sim.net.cells[gid].tags['cellType']
    netParams.stimSourceParams['IClamp_'+str(gid)] = { "type": "IClamp", "del": 20, "dur": sim_duration-20, "amp": 0.1 }
    netParams.stimSourceParams['EPSP_'+str(gid)] = { "type": "NetStim", "noise": 0.2, "start": 20, "rate": 40, "seed": int(seed) }
    netParams.stimSourceParams['IPSP_'+str(gid)] = { "type": "NetStim", "noise": 0.2, "start": 20, "rate": epsp_rate / eiratio, "seed": 20 }
    if stimulation_paradigm == 'iclampsoma':
    #this whole schema is built around the premise that there is one of each cell type whilch has its own stimultion
        netParams.stimTargetParams['IClamp->'+str(gid)] = {"conds":{"pop":cellPop,"cellType":cellType},"sec":sec_to_stim,"loc": 0.5,"source":'IClamp_'+str(gid)}
    elif stimulation_paradigm == 'epspsoma':
        netParams.stimTargetParams['EPSP->'+str(gid)] = {"source":"EPSP_"+str(gid),"conds":{"pop":cellPop,"cellType":cellType},"sec":sec_to_stim,"synMech":"AMPA","weight": 0.015}
    elif stimulation_paradigm == 'epspipspsoma':
        netParams.stimTargetParams['EPSP->'+str(gid)] = {"source":"EPSP_"+str(gid),"conds":{"pop":cellPop,"cellType":cellType},"sec":sec_to_stim,"synMech":"AMPA","weight": 0.015}
        netParams.stimTargetParams['IPSP->'+str(gid)] = {"source":"IPSP_"+str(gid),"conds":{"pop":cellPop,"cellType":cellType},"sec":sec_to_stim,"synMech":"GABA","weight": 0.015}
    elif stimulation_paradigm == 'epspipspsoma' or stimulation_paradigm == 'epspipspdendrites':
        for s,sec in enumerate(secs_to_stim_epsp):
            #add stims to every layer
            netParams.stimSourceParams['EPSP_'+str(gid)+'_'+sec] = { "type": "NetStim", "rate": epsp_rate/1000, "noise": 0.2, "start":1,"seed":s }
            netParams.stimTargetParams['EPSP->'+str(gid)+'_'+sec] = {"source":"EPSP_"+str(gid)+'_'+sec,"conds":{"pop":cellPop,"cellType":cellType},"sec":sec,"synMech":"AMPA","weight": 0.015,"delay":5}
        
        for s,sec in enumerate(secs_to_stim_ipsp):
            netParams.stimSourceParams['IPSP_'+str(gid)+'_'+sec] = { "type": "NetStim", "rate": epsp_rate / eiratio/1000, "noise": 0.2, "start": 1,"seed":s+1000 }
            netParams.stimTargetParams['IPSP->'+str(gid)+'_'+sec] = {"source":"IPSP_"+str(gid)+'_'+sec,"conds":{"pop":cellPop,"cellType":cellType},"sec":sec,"synMech":"GABA","weight": 0.015,"delay":5}
        

sim.create(netParams = netParams, simConfig = simConfig)

if stimulation_paradigm == 'epspipspdendrites':
    print("epsp stim:",secs_to_stim_epsp)
    print("ipsp stim:",secs_to_stim_ipsp)
    epspid = '_'.join([x.split('_')[1] for x in secs_to_stim_epsp])
    ipspid = '_'.join([x.split('_')[1] for x in secs_to_stim_ipsp])
    simid = '_'.join([epspid,ipspid])
    

#set variant values - only sets types listed below
for g,gid in enumerate(record_gids):
    if not sim.net.cells[gid].tags["cellType"] in ["INT-WT","INT-I1356M","INT-L479P","INT-KO"]:
        continue
    sectionlist = sim.net.cells[gid].secs.keys()
    cell_type = sim.net.cells[gid].tags["cellType"]

    dendlist = [s for s in sectionlist if "dend" in s]
    somalist = [s for s in sectionlist if "soma" in s]
    axonlist = [s for s in sectionlist if "axon" in s]

    secs_to_modify = dendlist.copy()
    secs_to_modify.extend(somalist)
    secs_to_modify.extend(axonlist)

    #changes in parameters for variant.  First row are parameter names, other rows are paraemters.  conductance and slopes (ms,hs) are multipliers, mh is an absolute shift
    data = {'name':['gNav11bar','mh','ms','hs','gNap_Et2bar'],
            'INT-WT':[1,0,1,1,1],
            'INT-I1356M':[0.56,-3.24,1.15,0.95,1],
            'INT-L479P':[1,0,1,1,0.6],
            'INT-KO':[0,0,1,1,1]}

    param_df = pd.DataFrame(data)
    for index, row in param_df.iterrows():
        for s,sec in enumerate(secs_to_modify):
            if row['name'] == 'gNap_Et2bar':
                paramVal = sim.net.cells[gid].secs[sec]['mechs']['Nap_Et2'][row['name']]
            elif row['name'] == 'gNav11bar':
                paramVal = sim.net.cells[gid].secs[sec]['mechs']['Nav11'][row['name']]
            else:
                paramVal = sim.net.cells[gid].secs[sec]['mechs']['Nav11var'][row['name']]

            if paramVal:
                if row['name'] == 'mh':
                    sim.net.modifyCells({'conds':{'cellType':cell_type},'secs':{sec:{'mechs':{'Nav11var':{row['name']:paramVal + row[cell_type]}}}}})
                elif row['name'] == 'gNav11bar':
                    #set the variant conductance
                    sim.net.modifyCells({'conds':{'cellType':cell_type},'secs':{sec:{'mechs':{'Nav11var':{'gNav11bar':paramVal * 0.5 * row[cell_type]}}}}})
                    #we also have to set the non-variant Nav11 conductance to 50%, the original code allowed looking at different percentages of each
                    sim.net.modifyCells({'conds':{'cellType':cell_type},'secs':{sec:{'mechs':{'Nav11':{'gNav11bar':paramVal * 0.5}}}}})
                elif row['name'] == 'gNap_Et2bar':
                    sim.net.modifyCells({'conds':{'cellType':cell_type},'secs':{sec:{'mechs':{'Nap_Et2':{row['name']:paramVal * row[cell_type]}}}}})
                else:
                    sim.net.modifyCells({'conds':{'cellType':cell_type},'secs':{sec:{'mechs':{'Nav11var':{row['name']:paramVal * row[cell_type]}}}}})


fig = plt.figure(figsize=(7.33,6.6),dpi=900) #15,9
subfigs = fig.subfigures(2,1,height_ratios=[2,1])
axs_left = subfigs[0].subplots(len(record_gids),1) 

results_dict = {new_list: [] for new_list in record_gids}
results_firstspike_dict = {}
results_firstdepoblock_dict = {}

#now set up everything for the run of simulations
#assuming the transition for the first spike is between 0 and initval
if stimulation_paradigm == 'iclampsoma':
    stopstepsize = 0.0001
    maxvaltotest = 10 
else:
    stopstepsize = .001  
    maxvaltotest = 10000  
    
initval = maxvaltotest 
maxval = {}
minval = {}
stepsize = {}
testval = {}
numberspikes = {}
lastnumberspikes = {}
maxfreq = {}
maxnumberspikes = {}
firstspike = {}
depofreq = {}
firstRun = {}
dealWithZeroStart = {}
for g,gid in enumerate(record_gids):
    maxval[g] = initval
    minval[g] = 0
    stepsize[g] = (maxval[g] - minval[g]) / 2

    testval[g] = initval
    lastnumberspikes[g] = 1
    maxfreq[g] = 0 #freq at which the number of spikes is the highest
    maxnumberspikes[g] = 1
    firstspike[g] = -1

    firstRun[g] = True
    dealWithZeroStart[g] = False

#first find the point at which spikes start... 
#here, we know that stepsizes for each neuron are the same, so we can just use the first...
while max(stepsize.values()) > stopstepsize:
    #find the next test value for each neuron
    for g,gid in enumerate(record_gids):
        if lastnumberspikes[g] >= maxnumberspikes[g]:
            maxnumberspikes[g] = lastnumberspikes[g]
            maxfreq[g] = testval[g]
        
        if lastnumberspikes[g] > 0 or dealWithZeroStart[g]:
            testval[g] = testval[g] - stepsize[g]
        else:
            testval[g] = testval[g] + stepsize[g]
        stepsize[g] /= 2
        testval[g] = np.round(testval[g],10)
        
    #set the next test value based on the stimulation paradigm
        if stimulation_paradigm == 'iclampsoma':
            sim.net.modifyStims({'conds':{'source':'IClamp_'+str(gid)},'amp':testval[g]})
        elif stimulation_paradigm == 'epspsoma':
            sim.net.modifyStims({'conds':{'source':'EPSP_'+str(gid)},'rate':testval[g]/1000})
        elif stimulation_paradigm == 'epspipspsoma':
            sim.net.modifyStims({'conds':{'source':'EPSP_'+str(gid)},'rate':testval[g]/1000})
            sim.net.modifyStims({'conds':{'source':'IPSP_'+str(gid)},'rate':testval[g]/eiratio/1000})
        elif stimulation_paradigm == 'epspipspsoma' or stimulation_paradigm == 'epspipspdendrites':
            for sec in secs_to_stim_epsp:
                sim.net.modifyStims({'conds':{'source':'EPSP_'+str(gid)+'_'+sec},'rate':testval[g]/1000})
            for sec in secs_to_stim_ipsp:
                sim.net.modifyStims({'conds':{'source':'IPSP_'+str(gid)+'_'+sec},'rate':testval[g]/eiratio/1000})

    sim.simulate()
    sim.analyze()

    data = {} 

    for g,gid in enumerate(record_gids):
        counts = np.unique(sim.allSimData['spkid'],return_counts=True)
        if np.any(counts) and gid in counts[0]:
            index = np.where(counts[0]==gid)
            #add tuple (testval, numspikes) for cell gid
            lastnumberspikes[g] = counts[1][index][0]
            if firstspike[g] == -1 or testval[g] < firstspike[g]:
                firstspike[g] = testval[g]
            results_dict[gid].append((testval[g],lastnumberspikes[g]))
        else:
            lastnumberspikes[g] = 0
            results_dict[gid].append((testval[g],0))

        data[g] = sim.allSimData['V_soma']['cell_'+str(gid)]
        for d in range(len(data[g])):
            data[g][d] = min(data[g][d],70)

    #deal with boundary case for very inhibited stimulation patterns
        if firstRun[g]:
            firstRun[g] = False
            if lastnumberspikes[g] == 0:
                dealWithZeroStart[g] = True
        if dealWithZeroStart[g]:
            if lastnumberspikes[g] > 0:
                dealWithZeroStart[g] = False

#now that we've found the first spike stimulus value, plot it for each neuron
for g,gid in enumerate(record_gids):
    times = np.linspace(0,sim.cfg.duration,len(data[g]))
    roundval = 2
    if stimulation_paradigm == 'iclampsoma':
        roundval = 4
    axs_left[g].plot(times,data[g],label='First spike at '+str(np.round(testval[g],roundval))+' Hz')
    if g==len(record_gids)-1:
        axs_left[g].set_xlabel("Time (ms)")
    else:
        plt.setp(axs_left[g].get_xticklabels(), visible=False)
    
    axs_left[g].set_yticks([-50,0])
    axs_left[g].title.set_text(sim.net.cells[gid].tags['pop'])
    axs_left[g].legend(loc='center left',bbox_to_anchor=(1.0,0.5))
    
    results_firstspike_dict[gid] = firstspike[g] 
    numberspikes[g] = 1 

    testval[g] = maxfreq[g]
    numberspikes[g] = maxnumberspikes[g]

#next, find a point beyond the depolarization block transition
while min(stepsize.values()) < maxvaltotest / 2:
    for g,gid in enumerate(record_gids):
        lastnumberspikes[g] = numberspikes[g]
        if numberspikes[g] >= maxnumberspikes[g]:
            maxnumberspikes[g] = numberspikes[g]
            maxfreq[g] = testval[g]

        stepsize[g] *= 2 
        testval[g] = testval[g] + stepsize[g]
 
    #set the next test value based on the stimulation paradigm
        if stimulation_paradigm == 'iclampsoma':
            sim.net.modifyStims({'conds':{'source':'IClamp_'+str(gid)},'amp':testval[g]})
        elif stimulation_paradigm == 'epspsoma':
            sim.net.modifyStims({'conds':{'source':'EPSP_'+str(gid)},'rate':testval[g]/1000})
        elif stimulation_paradigm == 'epspipspsoma':
            sim.net.modifyStims({'conds':{'source':'EPSP_'+str(gid)},'rate':testval[g]/1000})
            sim.net.modifyStims({'conds':{'source':'IPSP_'+str(gid)},'rate':testval[g]/eiratio/1000})
        elif stimulation_paradigm == 'epspipspsoma' or stimulation_paradigm == 'epspipspdendrites':
            for sec in secs_to_stim_epsp:
                sim.net.modifyStims({'conds':{'source':'EPSP_'+str(gid)+'_'+sec},'rate':testval[g]/1000})
            for sec in secs_to_stim_ipsp:
                sim.net.modifyStims({'conds':{'source':'IPSP_'+str(gid)+'_'+sec},'rate':testval[g]/eiratio/1000})

    sim.simulate()
    sim.analyze()

    data = {} 

    for g,gid in enumerate(record_gids):
        counts = np.unique(sim.allSimData['spkid'],return_counts=True)
        if np.any(counts) and gid in counts[0]:
            index = np.where(counts[0]==gid)
            #add tuple (testval, numspikes) for cell gid
            numberspikes[g] = counts[1][index][0]
            results_dict[gid].append((testval[g],numberspikes[g]))
        else:
            numberspikes[g] = 0
            results_dict[gid].append((testval[g],0))

        data[g] = sim.allSimData['V_soma']['cell_'+str(gid)]
        for d in range(len(data[g])):
            data[g][d] = min(data[g][d],70)

        depofreq[g] = testval[g]

#then do 10 fixed points 
testvals = np.array([1,5,10,30,50,100,150,200,250,300])
testvaluearray = np.zeros((len(record_gids),len(testvals)))
for g,gid in enumerate(record_gids):
    testvaluearray[g,:] = np.round(testvals,2)

#iterate over each column, which is a set of stimulation values, one for each neuron
for vals in testvaluearray.T:
    for g,gid in enumerate(record_gids):
        if stimulation_paradigm == 'iclampsoma':
            sim.net.modifyStims({'conds':{'source':'IClamp_'+str(gid)},'amp':vals[g]})
        elif stimulation_paradigm == 'epspsoma':
            sim.net.modifyStims({'conds':{'source':'EPSP_'+str(gid)},'rate':vals[g]/1000})
        elif stimulation_paradigm == 'epspipspsoma':
            sim.net.modifyStims({'conds':{'source':'EPSP_'+str(gid)},'rate':vals[g]/1000})
            sim.net.modifyStims({'conds':{'source':'IPSP_'+str(gid)},'rate':vals[g]/eiratio/1000})
        elif stimulation_paradigm == 'epspipspsoma' or stimulation_paradigm == 'epspipspdendrites':
            for sec in secs_to_stim_epsp:
                sim.net.modifyStims({'conds':{'source':'EPSP_'+str(gid)+'_'+sec},'rate':vals[g]/1000})
            for sec in secs_to_stim_ipsp:
                sim.net.modifyStims({'conds':{'source':'IPSP_'+str(gid)+'_'+sec},'rate':vals[g]/eiratio/1000})

    sim.simulate()
    sim.analyze()
        
    for g,gid in enumerate(record_gids):
        counts = np.unique(sim.allSimData['spkid'],return_counts=True)
        if np.any(counts) and gid in counts[0]:
            index = np.where(counts[0]==gid)
            #add tuple (testval, numspikes) for cell gid
            numberspikes[g] = counts[1][index][0]
            results_dict[gid].append((vals[g],numberspikes[g]))
        else:
            numberspikes[g] = 0
            results_dict[gid].append((vals[g],0))
 
        data[g] = sim.allSimData['V_soma']['cell_'+str(gid)]
        for d in range(len(data[g])):
            data[g][d] = min(data[g][d],70)
 
        if numberspikes[g] >= maxnumberspikes[g]:
            maxnumberspikes[g] = numberspikes[g]
            maxfreq[g] = vals[g]

#Finally, find the depolarization block transition (look between the point beyond (testval) and the max eifreq identified (maxfreq)
#here we have different stepsizes for each neuron...
for g,gid in enumerate(record_gids):
    maxval[g] = testval[g]
    minval[g] = maxfreq[g] 
    stepsize[g] = (maxval[g] - minval[g]) / 2

while max(stepsize.values()) > stopstepsize:
    for g,gid in enumerate(record_gids):
        if numberspikes[g] > maxnumberspikes[g]:
            maxnumberspikes[g] = numberspikes[g] #this is probably not good if we're here
            maxfreq[g]= testval[g]
            testval[g] = testval[g] + stepsize[g]
        else: 
            if testval[g] < maxfreq[g]:
                testval[g] = testval[g] + stepsize[g]
            else:
                testval[g] = testval[g] - stepsize[g]

        stepsize[g] /= 2
        testval[g] = np.round(testval[g],10)

        if stimulation_paradigm == 'iclampsoma':
            sim.net.modifyStims({'conds':{'source':'IClamp_'+str(gid)},'amp':testval[g]})
        elif stimulation_paradigm == 'epspsoma':
            sim.net.modifyStims({'conds':{'source':'EPSP_'+str(gid)},'rate':testval[g]/1000})
        elif stimulation_paradigm == 'epspipspsoma':
            sim.net.modifyStims({'conds':{'source':'EPSP_'+str(gid)},'rate':testval[g]/1000})
            sim.net.modifyStims({'conds':{'source':'IPSP_'+str(gid)},'rate':testval[g]/eiratio/1000})
        elif stimulation_paradigm == 'epspipspsoma' or stimulation_paradigm == 'epspipspdendrites':
            for sec in secs_to_stim_epsp:
                sim.net.modifyStims({'conds':{'source':'EPSP_'+str(gid)+'_'+sec},'rate':testval[g]/1000})
            for sec in secs_to_stim_ipsp:
                sim.net.modifyStims({'conds':{'source':'IPSP_'+str(gid)+'_'+sec},'rate':testval[g]/eiratio/1000})
        
    sim.simulate()
    sim.analyze()
 
    for g,gid in enumerate(record_gids):
        counts = np.unique(sim.allSimData['spkid'],return_counts=True)
        if np.any(counts) and gid in counts[0]:
            index = np.where(counts[0]==gid)
            #add tuple (testval, numspikes) for cell gid
            numberspikes[g] = counts[1][index][0]
            results_dict[gid].append((testval[g],numberspikes[g]))
        else:
            numberspikes[g] = 0
            results_dict[gid].append((testval[g],0))
 
        data[g] = sim.allSimData['V_soma']['cell_'+str(gid)]
        for d in range(len(data[g])):
            data[g][d] = min(data[g][d],70)



for g,gid in enumerate(record_gids):
    results_firstdepoblock_dict[gid] = testval[g]
     #plot max firing rate
    times = np.linspace(0,sim.cfg.duration,len(data[g]))
    if stimulation_paradigm == 'iclampsoma':
        axs_left[g].plot(times,data[g],label='Most spikes at '+str(np.round(testval[g],4))+' Hz')
    else:
        axs_left[g].plot(times,data[g],label='Most spikes at '+str(int(testval[g]))+' Hz')
    axs_left[g].set_yticks([-50,0])
    axs_left[g].title.set_text(sim.net.cells[gid].tags['pop'])
    axs_left[g].legend(loc='center left',bbox_to_anchor=(1.0,0.5))

#plot depo block
for g,gid in enumerate(record_gids):
    if stimulation_paradigm == 'iclampsoma':
        sim.net.modifyStims({'conds':{'source':'IClamp_'+str(gid)},'amp':depofreq[g]})
    elif stimulation_paradigm == 'epspsoma':
        sim.net.modifyStims({'conds':{'source':'EPSP_'+str(gid)},'rate':depofreq[g]/1000})
    elif stimulation_paradigm == 'epspipspsoma':
        sim.net.modifyStims({'conds':{'source':'EPSP_'+str(gid)},'rate':depofreq[g]/1000})
        sim.net.modifyStims({'conds':{'source':'IPSP_'+str(gid)},'rate':depofreq[g]/eiratio/1000})
    elif stimulation_paradigm == 'epspipspsoma' or stimulation_paradigm == 'epspipspdendrites':
        for sec in secs_to_stim_epsp:
            sim.net.modifyStims({'conds':{'source':'EPSP_'+str(gid)+'_'+sec},'rate':depofreq[g]/1000})
        for sec in secs_to_stim_ipsp:
            sim.net.modifyStims({'conds':{'source':'IPSP_'+str(gid)+'_'+sec},'rate':depofreq[g]/eiratio/1000})

    sim.simulate()
    sim.analyze()

    data[g] = sim.allSimData['V_soma']['cell_'+str(gid)]
    times = np.linspace(0,sim.cfg.duration,len(data[g]))
    if stimulation_paradigm == 'iclampsoma':
        axs_left[g].plot(times,data[g],label='Max EPSP tested: '+str(np.round(depofreq[g],4))+' Hz')
    else:
        axs_left[g].plot(times,data[g],label='Max EPSP tested: '+str(int(depofreq[g]))+' Hz')
    axs_left[g].set_yticks([-50,0])
    axs_left[g].legend(loc='center left',bbox_to_anchor=(1.0,0.5))


#write data to excel 
summary_df = pd.DataFrame(index=record_gids,columns=["cell type","first spike","first depolarization block","epsp_compartments","ipsp_compartments", "synapse stim arrangement type","synapse stim rating","synapse stim distance rating","dendrite stim categorization","randomization seed"])

#generate column names for the data 
column_names = []
for i in record_gids:
    cell = sim.net.cells[i].tags['cellType']
    column_names.append(cell + " freq")
    column_names.append(cell + " data")

data_df = pd.DataFrame(columns=column_names)

summary_df.at[0,"epsp_compartments"] = epspid
summary_df.at[0,"ipsp_compartments"] = ipspid
summary_df.at[0,"synapse stim rating"] = syn_exc_rating 
summary_df.at[0,"synapse stim distance rating"] = syn_dist_rating 
summary_df.at[0,"dendrite stim categorization"] = top_sec_exc_list 
summary_df.at[0,"randomization seed"] = seed 

outdir = ""

if syn_exc_rating < 0.0:
    summary_df.at[0,"synapse stim arrangement type"] = "inhibitory"
    outdir = "results_inh/"
elif syn_exc_rating > 2.0:
    summary_df.at[0,"synapse stim arrangement type"] = "excitatory"
    outdir = "results_exc/"
else:
    summary_df.at[0,"synapse stim arrangement type"] = "mixed"
    outdir = "results_mix/"

if putInTopDirectory:
    outdir = ""

for i in record_gids:
    summary_df.at[i,"cell type"] = sim.net.cells[i].tags['cellType']
    summary_df.at[i,"first spike"] = results_firstspike_dict[i]
    summary_df.at[i,"first depolarization block"] = results_firstdepoblock_dict[i]

    #for tup in results_dict[i]:
    cell = sim.net.cells[i].tags['cellType']
    freqlist = [m for m, n in results_dict[i]]
    datalist = [n for m, n in results_dict[i]]
    data_df[cell+" freq"] = pd.Series(freqlist)
    data_df[cell+" data"] = pd.Series(datalist)

if stimulation_paradigm == 'iclampsoma':
    outputxlsxfile = 'output_data_iclamprun.xlsx'
    outputpngfile = 'vary_iclamp.png'
elif stimulation_paradigm == 'epspsoma':
    outputxlsxfile = 'output_data_epsprun'+str(seed)+'.xlsx'
    outputpngfile = 'vary_epsp'+str(seed)+'.png'
elif stimulation_paradigm == 'epspipspsoma':
    outputxlsxfile = 'output_data_ei'+str(eiratio)+'_'+str(seed)+'.xlsx'
    outputpngfile = 'vary_epsp_ei' + str(eiratio)+'_'+str(seed)+'.png'
else:
    outputxlsxfile = 'output_data_ei'+str(eiratio)+'_id_' + simid + '.xlsx'
    outputpngfile = 'vary_EPSP_ei' + str(eiratio)+'id_' + simid + '.png'
with pd.ExcelWriter(outdir + outputxlsxfile) as writer:  
    summary_df.to_excel(writer, sheet_name='Summary')
    data_df.to_excel(writer, sheet_name='Data')

maxdepoblock = max(results_firstdepoblock_dict.values())
subfigs[0].text(0.003, 0.5, "Soma membrane potential (mV)", va='center',rotation='vertical')
subfigs[0].text(0.03, 0.96, "A)",fontsize=15)
subfigs[0].text(0.03, 0.73, "B)",fontsize=15)
subfigs[0].text(0.03, 0.51, "C)",fontsize=15)
subfigs[0].text(0.03, 0.29, "D)",fontsize=15)
subfigs[1].text(0.03, 0.98, "E)",fontsize=15)


lastax = len(record_gids)

#plot the summary figure
numspikesplot = subfigs[1].subplots(1,1)

for i in record_gids:
    results_dict[i].sort(key = lambda x: x[0])
    p = numspikesplot.plot(*zip(*results_dict[i]),marker='o',markersize = 4,alpha = 0.75,label=sim.net.cells[i].tags['pop'])
    if results_firstspike_dict[i] >= 0:
        numspikesplot.axvline(results_firstspike_dict[i],alpha = 0.75,color = p[-1].get_color())
    numspikesplot.axvline(results_firstdepoblock_dict[i],alpha = 0.75,color = p[-1].get_color())

numspikesplot.set_xlabel("EPSP Freqency (Hz)")
numspikesplot.set_ylabel("Number spikes")
if stimulation_paradigm == 'iclampsoma':
    numspikesplot.set_xlim([0,1])
else:
    numspikesplot.set_xlim([0,min(maxvaltotest,max(600,maxdepoblock * 1.5))])
numspikesplot.legend(loc='upper right')
numspikesplot.set_title("Genetic Variant Effect on Neuronal Firing",fontweight='bold')

subfigs[0].subplots_adjust(left=.075,right=.6,top=0.95, bottom=0.15,wspace=0.2,hspace=0.4)
subfigs[1].subplots_adjust(left=.075,right=.95,top=0.95,bottom=0.2,wspace=0.1,hspace=0.2)

plt.savefig(outdir + outputpngfile)

if stimulation_paradigm == 'epspipspdendrites':
    print("epsp stim:",secs_to_stim_epsp)
    print("ipsp stim:",secs_to_stim_ipsp)