HB-PAC_disinhibitory_network/simSynCurrents.py

import os
import sys

import seaborn as sns
from brian2 import defaultclock, units, StateMonitor
from matplotlib import pyplot as plt
from brian2.core.network import Network
from dirDefs import homeFolder
from models.neuronModels import VSNeuron, JOSpikes265, getSineInput
from models.neurons import AdExp
from models.synapses import exp2Syn, exp2SynStateInits
from mplPars import mplPars
from paramLists import synapsePropsList, inputParsList, AdExpPars


from neo import AnalogSignal
import nixio
from neoNIXIO import addAnalogSignal2Block, addMultiTag
import quantities as qu
from brianUtils import addBrianQuantity2Section
sns.set(style="whitegrid", rc=mplPars)


simSettleTime = 600 * units.ms

simStepSize = 0.1 * units.ms
simDuration = 150 * units.ms
# inputParsName = 'onePulse'
# inputParsName = 'twoPulse'
# inputParsName = 'threePulse'
inputParsName = "thirtyMSPulse"
# inputParsName = "fortyMSPulse"
showBefore = 50 * units.ms
showAfter = 50 * units.ms

# simStepSize = 0.1 * units.ms
# simDuration = 450 * units.ms
# # inputParsName = "pTShortInt20Dur10"
# # inputParsName = "pTShortInt20Dur16"
# # inputParsName = "pTShortInt33Dur10"
# # inputParsName = "pTShortInt33Dur16"
# # inputParsName = "pTShortInt33Dur20"
# # inputParsName = "pTShortInt50Dur10"
# # inputParsName = "pTShortInt50Dur16"
# # inputParsName = "pTShortInt50Dur20"
# inputParsName = "pTShortInt100Dur10"
# # inputParsName = "pTShortInt100Dur16"
# # inputParsName = "pTShortInt100Dur20"
#
# showBefore = 100 * units.ms
# showAfter = 100 * units.ms

# simStepSize = 0.1 * units.ms
# simDuration = 1500 * units.ms
# # inputParsName = 'oneSecondPulse'
# # inputParsName = 'pulseTrainInt20Dur10'
# inputParsName = 'pulseTrainInt20Dur16'
# # inputParsName = 'pulseTrainInt33Dur10'
# # inputParsName = 'pulseTrainInt33Dur16'
# showBefore = 500 * units.ms
# showAfter = 500 * units.ms

DLInt1ModelProps = "DLInt1Aynur"

DLInt1SynapsePropsE = 'DLInt1_syn_try2_e'
# DLInt1SynapsePropsE = ""
DLInt1SynapsePropsI = 'DLInt1_syn_try2_i'
# DLInt1SynapsePropsI = ""
DLInt1SynapseProps = "".join((DLInt1SynapsePropsE, DLInt1SynapsePropsI))

DLInt2ModelProps = "DLInt2Try2"

DLInt2SynapseProps = 'DLInt2_syn_try2'

DLInt1DLInt2SynProps = "DLInt1_DLInt2_try1"

opDir = os.path.join(homeFolder, DLInt1ModelProps + DLInt2ModelProps,
                         DLInt1SynapseProps + DLInt2SynapseProps + DLInt1DLInt2SynProps,
                         inputParsName)
opFileDLInt1 = os.path.join(opDir, 'SynCurrentTracesDLInt1.png')
opFileDLInt2 = os.path.join(opDir, 'SynCurrentTracesDLInt2.png')
OPNixFile = os.path.join(opDir, 'simResWithSynCurrents.h5')
if os.path.isfile(OPNixFile):
    ch = input('Results already exist at {}. Delete?(y/n):'.format(OPNixFile))
    if ch == 'y':
        os.remove(OPNixFile)
        if os.path.isfile(opFileDLInt1):
            os.remove(opFileDLInt1)
        if os.path.isfile(opFileDLInt2):
            os.remove(opFileDLInt2)
    else:
        sys.exit('User Abort!')

elif not os.path.isdir(opDir):
    os.makedirs(opDir)

inputPars = getattr(inputParsList, inputParsName)


net = Network()
JO = JOSpikes265(nOutputs=1, simSettleTime=simSettleTime, **inputPars)
net.add(JO.JOSGG)

DLInt1PropsDict = getattr(AdExpPars, DLInt1ModelProps)
dlint1 = VSNeuron(**AdExp, inits=DLInt1PropsDict, name='dlint1')
dlint1.recordSpikes()
dlint1.recordMembraneV()

if DLInt1SynapsePropsE:
    synPropsEDLInt1 = getattr(synapsePropsList, DLInt1SynapsePropsE)
    dlint1.addSynapse(synName="ExiJO", sourceNG=JO.JOSGG, **exp2Syn,
                          synParsInits=synPropsEDLInt1,
                          synStateInits=exp2SynStateInits,
                          sourceInd=0, destInd=0
                          )
if DLInt1SynapsePropsI:
    synPropsIDLInt1 = getattr(synapsePropsList, DLInt1SynapsePropsI)
    dlint1.addSynapse(synName="InhJO", sourceNG=JO.JOSGG, **exp2Syn,
                          synParsInits=synPropsIDLInt1,
                          synStateInits=exp2SynStateInits,
                          sourceInd=0, destInd=0
                          )

dlint1.addToNetwork(net)

if DLInt1SynapsePropsE:
    gEMonitorDLInt1 = StateMonitor(dlint1.incomingSynapses["ExiJO"], "g_ExiJO", record=True)
    net.add(gEMonitorDLInt1)

if DLInt1SynapsePropsI:
    gIMonitorDLInt1 = StateMonitor(dlint1.incomingSynapses["InhJO"], "g_InhJO", record=True)
    net.add(gIMonitorDLInt1)

DLInt2PropsDict = getattr(AdExpPars, DLInt2ModelProps)
dlint2 = VSNeuron(**AdExp, inits=DLInt2PropsDict, name='dlint2')
dlint2.recordMembraneV()
dlint2.recordSpikes()

if DLInt2SynapseProps:
    synPropsEDLInt2 = getattr(synapsePropsList, DLInt2SynapseProps)
    dlint2.addSynapse(synName="ExiJO", sourceNG=JO.JOSGG, **exp2Syn,
                          synParsInits=synPropsEDLInt2,
                          synStateInits=exp2SynStateInits,
                          sourceInd=0, destInd=0
                          )

if DLInt1DLInt2SynProps:
    synPropsIDLInt2 = getattr(synapsePropsList, DLInt1DLInt2SynProps)
    dlint2.addSynapse(synName="DLInt1", sourceNG=dlint1.ng, **exp2Syn,
                          synParsInits=synPropsIDLInt2,
                          synStateInits=exp2SynStateInits,
                          sourceInd=0, destInd=0
                          )


dlint2.addToNetwork(net)

if DLInt2SynapseProps:
    gEMonitorDLInt2 = StateMonitor(dlint2.incomingSynapses["ExiJO"], "g_ExiJO", record=True)
    net.add(gEMonitorDLInt2)

if DLInt1DLInt2SynProps:
    gIMonitorDLInt2 = StateMonitor(dlint2.incomingSynapses["DLInt1"], "g_DLInt1", record=True)
    net.add(gIMonitorDLInt2)

defaultclock.dt = simStepSize
totalSimDur = simDuration + simSettleTime
net.run(totalSimDur, report='text')

simT, DLInt2memV = dlint2.getMemVTrace()
DLInt2spikeTimes = dlint2.getSpikes()

dlint2MemVAS = AnalogSignal(signal=DLInt2memV / units.mV,
                            sampling_period=(simStepSize / units.ms) * qu.ms,
                            t_start=0 * qu.mV,
                            units="mV",
                            name="DLInt2 MemV")

dlint2SpikesQU = (DLInt2spikeTimes / units.ms) * qu.ms

simT, DLInt1memV = dlint1.getMemVTrace()
DLInt1spikeTimes = dlint1.getSpikes()

dlint1MemVAS = AnalogSignal(signal=DLInt1memV / units.mV,
                            sampling_period=(simStepSize / units.ms) * qu.ms,
                            t_start=0 * qu.mV,
                            units="mV",
                            name="DLInt1 MemV")

dlint1SpikesQU = (DLInt1spikeTimes / units.ms) * qu.ms


joSpikesQU = (JO.spikeTimes / units.ms) * qu.ms

sineInput = getSineInput(simDur=simDuration, simStepSize=simStepSize,
                         sinPulseDurs=inputPars['sinPulseDurs'],
                         sinPulseStarts=inputPars['sinPulseStarts'],
                         freq=265 * units.Hz, simSettleTime=simSettleTime)

inputAS = AnalogSignal(signal=sineInput,
                                sampling_period=(simStepSize / units.ms) * qu.ms,
                                t_start=0 * qu.mV,
                                units="um",
                                name="Input Vibration Signal")

fig, axs = plt.subplots(nrows=3, figsize=(10, 6.25), sharex='col')
axs[0].plot(simT / units.ms, DLInt2memV / units.mV)
spikesY = DLInt2memV.min() + 1.05 * (DLInt2memV.max() - DLInt2memV.min())
axs[0].plot(DLInt2spikeTimes / units.ms, [spikesY / units.mV] * DLInt2spikeTimes.shape[0], 'k^')
axs[0].set_ylabel('DLInt2\nMembrane\nPotential\n(mV)')
axs[0].set_xlim([(simSettleTime - showBefore) / units.ms,
                     (totalSimDur + showAfter) / units.ms])

fig1, axs1 = plt.subplots(nrows=3, figsize=(10, 6.25), sharex='col')
axs1[0].plot(simT / units.ms, DLInt1memV / units.mV)
spikesY = DLInt1memV.min() + 1.05 * (DLInt1memV.max() - DLInt1memV.min())
axs1[0].plot(DLInt1spikeTimes / units.ms, [spikesY / units.mV] * DLInt1spikeTimes.shape[0], 'k^')
axs1[0].set_ylabel('DLInt1\nMembrane\nPotential\n(mV)')
axs1[0].set_xlim([(simSettleTime - showBefore) / units.ms,
                     (totalSimDur + showAfter) / units.ms])

if DLInt2SynapseProps:
    gSynEDLInt2 = gEMonitorDLInt2.g_ExiJO[0]
    iSynEDLInt2 = -gSynEDLInt2 * (DLInt2memV - synPropsEDLInt2['Esyn'])
    axs[1].plot(simT / units.ms,
                iSynEDLInt2 / units.nA, 'm-', label=r'$I_{synE}$')
    iSynEASDLInt2 = AnalogSignal(signal=iSynEDLInt2 / units.nA,
                           sampling_period=(simStepSize / units.ms) * qu.ms,
                           t_start=0 * qu.mV,
                           units="nA",
                           name="DL-Int-2 input EPSC")

if DLInt1DLInt2SynProps:
    gSynI = gIMonitorDLInt2.g_DLInt1[0]
    iSynI = -gSynI * (DLInt2memV - synPropsIDLInt2['Esyn'])
    axs[1].plot(simT / units.ms,
                iSynI / units.nA, 'c-', label=r'$I_{synI}$')
    iSynIASDLInt2 = AnalogSignal(signal=iSynI / units.nA,
                           sampling_period=(simStepSize / units.ms) * qu.ms,
                           t_start=0 * qu.mV,
                           units="nA",
                           name="DL-Int-2 input IPSC")

axs[1].legend(loc='center right')
axs[1].set_ylabel("Synaptic\ncurrents\n(nA)")

axs[2].plot(simT / units.ms, sineInput,
            color=[130. / 255, 72. / 255, 7. / 255], ls='-', marker='None',
            label='Vibration Input')
axs[2].plot(JO.spikeTimes / units.ms, [sineInput.max() * 1.05] * len(JO.spikeTimes), 'k^',
        label='JO Spikes')
axs[2].legend(loc='center right')
axs[2].set_xlabel('time (ms)')
axs[2].set_ylabel('Input')
fig.tight_layout()
fig.canvas.draw()
# plt.show()
fig.savefig(opFileDLInt2, dpi=150)

if DLInt1SynapsePropsE:
    gSynEDLInt1 = gEMonitorDLInt1.g_ExiJO[0]
    iSynEDLInt1 = -gSynEDLInt1 * (DLInt1memV - synPropsEDLInt1['Esyn'])
    axs1[1].plot(simT / units.ms,
                iSynEDLInt1 / units.nA, 'r-', label=r'$I_{synE}$')
    iSynEASDLInt1 = AnalogSignal(signal=iSynEDLInt1 / units.nA,
                                sampling_period=(simStepSize / units.ms) * qu.ms,
                                t_start=0 * qu.mV,
                                units="nA",
                                name="DL-Int-1 input EPSC")

if DLInt1SynapsePropsI:
    gSynIDLInt1 = gIMonitorDLInt1.g_InhJO[0]
    iSynIDLInt1 = -gSynIDLInt1 * (DLInt1memV - synPropsIDLInt1['Esyn'])
    axs1[1].plot(simT / units.ms,
                iSynIDLInt1 / units.nA, 'g-', label=r'$I_{synI}$')
    iSynIASDLInt1 = AnalogSignal(signal=iSynIDLInt1 / units.nA,
                           sampling_period=(simStepSize / units.ms) * qu.ms,
                           t_start=0 * qu.mV,
                           units="nA",
                           name="DL-Int-1 input IPSC")


axs1[1].legend(loc='center right')
axs1[1].set_ylabel("Synaptic\ncurrents\n(nA)")

axs1[2].plot(simT / units.ms, sineInput,
            color=[130. / 255, 72. / 255, 7. / 255], ls='-', marker='None',
            label='Vibration Input')
axs1[2].plot(JO.spikeTimes / units.ms, [sineInput.max() * 1.05] * len(JO.spikeTimes), 'k^',
        label='JO Spikes')
axs1[2].legend(loc='center right')
axs1[2].set_xlabel('time (ms)')
axs1[2].set_ylabel('Input')
fig1.tight_layout()
fig1.canvas.draw()
# plt.show()
fig1.savefig(opFileDLInt1, dpi=150)

nixFile = nixio.File.open(OPNixFile, mode=nixio.FileMode.ReadWrite)
neuronModels = nixFile.create_section("Neuron Models", "Model Parameters")

DLInt1PropsSec = neuronModels.create_section("DL-Int-1", "AdExp")

for propName, propVal in DLInt1PropsDict.items():
    addBrianQuantity2Section(DLInt1PropsSec, propName, propVal)

DLInt2PropsSec = neuronModels.create_section("DL-Int-2", "AdExp")

for propName, propVal in DLInt2PropsDict.items():
    addBrianQuantity2Section(DLInt2PropsSec, propName, propVal)

inputSec = nixFile.create_section("Input Parameters", "Sinusoidal Pulses")

for parName, parVal in inputPars.items():
    addBrianQuantity2Section(inputSec, parName, parVal)

addBrianQuantity2Section(inputSec, "simSettleTime", simSettleTime)

brianSimSettingsSec = nixFile.create_section("Simulation Parameters", "Brian Simulation")
addBrianQuantity2Section(brianSimSettingsSec, "simStepSize", simStepSize)
addBrianQuantity2Section(brianSimSettingsSec, "totalSimDuration", totalSimDur)
brianSimSettingsSec.create_property("method", nixio.Value("euler"))


synPropsSec = nixFile.create_section("Synapse Models", "Model Parameters")

if DLInt1SynapsePropsE:

    JODLInt1SynESec = synPropsSec.create_section("JODLInt1Exi", "DoubleExpSyn")
    JODLInt1SynEDict = getattr(synapsePropsList, DLInt1SynapsePropsE)

    for propName, propVal in JODLInt1SynEDict.items():
        addBrianQuantity2Section(JODLInt1SynESec, propName, propVal)

    JODLInt1SynESec.create_property("PreSynaptic Neuron", nixio.Value("JO"))
    JODLInt1SynESec.create_property("PostSynaptic Neuron", nixio.Value("DLInt1"))


if DLInt1SynapsePropsI:

    JODLInt1SynISec = synPropsSec.create_section("JODLInt1Inh", "DoubleExpSyn")
    JODLInt1SynIDict = getattr(synapsePropsList, DLInt1SynapsePropsI)

    for propName, propVal in JODLInt1SynIDict.items():
        addBrianQuantity2Section(JODLInt1SynISec, propName, propVal)
    JODLInt1SynISec.create_property("PreSynaptic Neuron", nixio.Value("JO"))
    JODLInt1SynISec.create_property("PostSynaptic Neuron", nixio.Value("DLInt1"))

if DLInt2SynapseProps:

    JODLInt2SynESec = synPropsSec.create_section("JODLInt2Exi", "DoubleExpSyn")
    JODLInt2SynEDict = getattr(synapsePropsList, DLInt2SynapseProps)

    for propName, propVal in JODLInt2SynEDict.items():
        addBrianQuantity2Section(JODLInt2SynESec, propName, propVal)
    JODLInt2SynESec.create_property("PreSynaptic Neuron", nixio.Value("JO"))
    JODLInt2SynESec.create_property("PostSynaptic Neuron", nixio.Value("DLInt2"))

if DLInt1DLInt2SynProps:

    DLInt1DLInt2SynSec = synPropsSec.create_section("DLInt1DLInt2Inh", "DoubleExpSyn")
    DLInt1DLInt2SynDict = getattr(synapsePropsList, DLInt1DLInt2SynProps)

    for propName, propVal in DLInt1DLInt2SynDict.items():
        addBrianQuantity2Section(DLInt1DLInt2SynSec, propName, propVal)
    DLInt1DLInt2SynSec.create_property("PreSynaptic Neuron", nixio.Value("DLInt1"))
    DLInt1DLInt2SynSec.create_property("PostSynaptic Neuron", nixio.Value("DLInt2"))

blk = nixFile.create_block("Simulation Traces", "Brian Output")
DLInt2DA = addAnalogSignal2Block(blk, dlint2MemVAS)
DLInt1DA = addAnalogSignal2Block(blk, dlint1MemVAS)
inputDA = addAnalogSignal2Block(blk, inputAS)
if DLInt1SynapsePropsE:
    epscAS = addAnalogSignal2Block(blk, iSynEASDLInt1)
if DLInt1SynapsePropsI:
    ipscAS = addAnalogSignal2Block(blk, iSynIASDLInt1)
if DLInt2SynapseProps:
    epscAS = addAnalogSignal2Block(blk, iSynEASDLInt2)
if DLInt1DLInt2SynProps:
    ipscAS = addAnalogSignal2Block(blk, iSynIASDLInt2)

addMultiTag("DLInt2 Spikes", type="Spikes", positions=dlint2SpikesQU,
            blk=blk, refs=[DLInt2DA])
addMultiTag("DLInt1 Spikes", type="Spikes", positions=dlint1SpikesQU,
            blk=blk, refs=[DLInt2DA])
addMultiTag("JO Spikes", type="Spikes", positions=joSpikesQU,
            blk=blk, refs=[inputDA])

nixFile.close()