<reponame>subhacom/mbnet
# kc_ggn_feedback_dclamp.py --- 
# Author: <NAME>
# Created: Tue Aug 20 10:58:08 2019 (-0400)
# Last-Updated: Wed Dec 11 17:32:49 2019 (-0500)
#           By: <NAME>
# Version: $Id$

# Code:
"""This script for testing expansion of the dynamic range of a KC due to GGN inhibition.

Instead of running a whole simulation in the full network, we play the
GGN membrane potential back to the KC.

We use GGN Vm from two simulations, one with low PN activity and
another with high PN activity.

"""
from __future__ import print_function
import os
import sys
sys.path += ['D:/subhasis_ggn/model/mb', 'D:/subhasis_ggn/model/mb/network', 'D:/subhasis_ggn/model/nrn']
import argparse
import numpy as np
from collections import defaultdict
import matplotlib
matplotlib.use('Qt5Agg')
from matplotlib import pyplot as plt
import h5py as h5
from config import Q_, h, logger, timestamp, mypid, myjobid, nrn_version
from timeit import default_timer as timer
import ephys
import nrnutils as nu
import neurograph as ng
import nsdf


GGN_KC_SYN_PARAMS = {
    'vmid': Q_('-40mV').to('mV').m,
    'vslope': Q_('5.0mV').to('mV').m,
    'e': Q_('-80mV').to('mV').m,
    'gbar': Q_('1e-3uS').to('uS').m,
    'tau': Q_('4.0ms').to('ms').m
}

# keep global reference of created model components so that they are
# not garbage collected when out of scope
model_dict = {}

def make_kc_with_dynaclamp(kc_name, kc_file, inject, tstart, tend, ggn_vm=None):
    """Read KC model from `kc_file`, inject current `inject` nA, apply
    dynamic clamp `ggn_vm`, which should be a 2D array with time (ms)
    in column 0, and voltage (mV) in column 1.

    """
    global model_dict
    kc = nu.create_cell(kc_name, filename=kc_file)
    model_dict[kc] = None
    iclamp = ephys.setup_current_clamp(kc.soma, pos=0.5, delay=Q_(tstart, 'ms'),
                                       duration=Q_((tend - tstart), 'ms'),
                                       amplitude=Q_(inject, 'nA'))
    model_dict[iclamp] = None
    ggn_g_vec = None
    if ggn_vm is not None:
        syn = h.GradedSyn(kc.soma(0.5))
        for attr, value in GGN_KC_SYN_PARAMS.items():
            setattr(syn, attr, value)
        model_dict[syn] = None
        ggn_comp = h.Section('ggn')
        model_dict[ggn_comp] = None
        h.setpointer(ggn_comp(0.5)._ref_v, 'vpre', syn)
        ggn_vm_vec = h.Vector(ggn_vm[:, 1])
        tvec = h.Vector(ggn_vm[:, 0])
        model_dict[tvec] = None
        # vec.play(var_reference, t, continuous) for interpolating        
        ret = ggn_vm_vec.play(ggn_comp(0.5)._ref_v, tvec, 1)
        print('####', ret)
        model_dict[ggn_vm_vec] = None
        ggn_g_vec = h.Vector()
        ggn_g_vec.record(syn._ref_g)
        model_dict[ggn_g_vec] = None
    kc_vm_vec = h.Vector()
    kc_vm_vec.record(kc.soma(0.5)._ref_v)
    model_dict[kc_vm_vec] = None
    print('Built model')
    return (kc_vm_vec, ggn_g_vec)

def make_parser():
    parser = argparse.ArgumentParser(description='Simulate KC with GGN inhibition at multiple current injections')
    parser.add_argument('--kc-file', type=str, dest='kc_file',
                        required=True, help='KC cell template file'
                        ' (.hoc)')
    parser.add_argument('--kc', type=str, dest='kc', required=True,
                        help='KC cell template name in template file')
    parser.add_argument('--ggn-vm-file', type=str, action='append', dest='ggn_vm_file',
                        required=True,
                        help='CSV file with column 0 time in ms, column 1 GGN Vm in mV')
    parser.add_argument('--istart', type=str, dest='istart', help='Starting amplitude of current (with unit)')
    parser.add_argument('--iend', type=str, dest='iend', help='Ending amplitude of current (with unit)')
    parser.add_argument('--di', type=str, dest='di', help='Current increments (with unit)')
    parser.add_argument('--tstart', type=str, dest='tstart', help='Current injection start time (with unit)')
    parser.add_argument('--tend', type=str, dest='tend', help='Current injection end time (with unit)')    
    return parser


def main():
    parser = make_parser()
    args = parser.parse_args()
    logger.info('Command line args: {}'.format(str(sys.argv)))
    print(args.ggn_vm_file)
    # KCs with GGN inhibition
    inhibited_vec = defaultdict(list)
    solo_vec_list = []
    tstart = Q_(args.tstart).to('ms').m
    tend = Q_(args.tend).to('ms').m
    istart = Q_(args.istart).to('nA').m
    iend = Q_(args.iend).to('nA').m
    di = Q_(args.di).to('nA').m
    irange = np.arange(istart, iend + di/2.0, di)
    logger.info('Starting current: {} nA'.format(istart))
    logger.info('End current: {} nA'.format(iend))
    logger.info('Increment: {} nA'.format(di))
    logger.info('current range: {}'.format(irange))
    ggn_vm = {}
    for input_file in args.ggn_vm_file:
        ggn_vm[input_file] = np.loadtxt(input_file)
    for inject in irange:
        for input_file, vm in ggn_vm.items():
            kc_vvec, ggn_gvec = make_kc_with_dynaclamp(args.kc, args.kc_file, inject, tstart, tend, vm)
            inhibited_vec[input_file].append((kc_vvec, ggn_gvec))
        # KC without any inhibition
        kc_vvec, ggn_gvec = make_kc_with_dynaclamp(args.kc, args.kc_file, inject, tstart, tend)
        solo_vec_list.append(kc_vvec)
    tvec = h.Vector()
    tvec.record(h._ref_t)
    h.tstop = tend
    print('Init')    
    h.init()
    print('Run')
    h.run()
    print('Finished simulation')
    fig, ax = plt.subplots(nrows=len(irange)+1, ncols=len(ggn_vm)+1, sharex='all', sharey='all')
    t = np.array(tvec.x)
    solo_data = []
    for ii, vvec in enumerate(solo_vec_list):
        ax[ii+1, 0].plot(tvec, vvec, color='#e66101')
        solo_data.append(np.array(vvec.x))
    combined = np.vstack(solo_data)
    
    prefix = 'UTC' + timestamp.strftime('%Y%m%d_%H%M%S')
    fname = '{}_solo_kc.npz'.format(prefix)
    np.savez(fname,
             t=t,
             vm=combined,
             inject=irange)
    logger.info('Saved solo KC data in {}'.format(fname))
    for jj, input_file in enumerate(args.ggn_vm_file):
        fname = '{}_{}.npz'.format(prefix, os.path.basename(input_file))
        data = []
        kc_vm_list = inhibited_vec[input_file]
        for ii, (vvec, gvec) in enumerate(kc_vm_list):
            data.append(np.array(vvec.x))
            ax[ii+1, jj+1].plot(tvec, vvec, color='#e66101')
            ax[ii+1, 0].set_ylabel('{} pA'.format(irange[ii]*1e3))
            # ax[ii+1, 0].set_ylabel('{} pA'.format(int(np.round(irange[ii]*1e3))))  # to avoid decimal point when integer values 
        # ax[0, jj+1].plot(tvec, gvec)
        # ax[0, jj+1].plot(ggn_vm[input_file][:,0], ggn_vm[input_file][:,1])
        ax[0, jj+1].set_title(input_file)
        combined = np.vstack(data)
        np.savez(fname, combined=combined, irange=irange, ggn_vm=ggn_vm[input_file])
        logger.info('Saved data from dynamic clamp with input from {} in {}'.format(
            input_file, fname))
    for axis in ax.flat:
        axis.set_xlim(250, 1750)
        fig.set_size_inches(210/25.4, 290/25.4)
    fig.tight_layout()
    fig.savefig('{}_KC_dynamic_range_with_ggn_vm.svg'.format(prefix))
    plt.show()
    print('End')

    
if __name__ == '__main__':
    main()
    

# 
# kc_ggn_feedback_dclamp.py ends here

# kc_ggn_feedback_dclamp.py --- 
# Author: <NAME>
# Created: Tue Aug 20 10:58:08 2019 (-0400)
# Last-Updated: Wed Dec 11 17:32:49 2019 (-0500)
#           By: <NAME>
# Version: $Id$

# Code:
"""This script for testing expansion of the dynamic range of a KC due to GGN inhibition.

Instead of running a whole simulation in the full network, we play the
GGN membrane potential back to the KC.

We use GGN Vm from two simulations, one with low PN activity and
another with high PN activity.

"""
from __future__ import print_function
import os
import sys
sys.path += ['D:/subhasis_ggn/model/mb', 'D:/subhasis_ggn/model/mb/network', 'D:/subhasis_ggn/model/nrn']
import argparse
import numpy as np
from collections import defaultdict
import matplotlib
matplotlib.use('Qt5Agg')
from matplotlib import pyplot as plt
import h5py as h5
from config import Q_, h, logger, timestamp, mypid, myjobid, nrn_version
from timeit import default_timer as timer
import ephys
import nrnutils as nu
import neurograph as ng
import nsdf


GGN_KC_SYN_PARAMS = {
    'vmid': Q_('-40mV').to('mV').m,
    'vslope': Q_('5.0mV').to('mV').m,
    'e': Q_('-80mV').to('mV').m,
    'gbar': Q_('1e-3uS').to('uS').m,
    'tau': Q_('4.0ms').to('ms').m
}

# keep global reference of created model components so that they are
# not garbage collected when out of scope
model_dict = {}

def make_kc_with_dynaclamp(kc_name, kc_file, inject, tstart, tend, ggn_vm=None):
    """Read KC model from `kc_file`, inject current `inject` nA, apply
    dynamic clamp `ggn_vm`, which should be a 2D array with time (ms)
    in column 0, and voltage (mV) in column 1.

    """
    global model_dict
    kc = nu.create_cell(kc_name, filename=kc_file)
    model_dict[kc] = None
    iclamp = ephys.setup_current_clamp(kc.soma, pos=0.5, delay=Q_(tstart, 'ms'),
                                       duration=Q_((tend - tstart), 'ms'),
                                       amplitude=Q_(inject, 'nA'))
    model_dict[iclamp] = None
    ggn_g_vec = None
    if ggn_vm is not None:
        syn = h.GradedSyn(kc.soma(0.5))
        for attr, value in GGN_KC_SYN_PARAMS.items():
            setattr(syn, attr, value)
        model_dict[syn] = None
        ggn_comp = h.Section('ggn')
        model_dict[ggn_comp] = None
        h.setpointer(ggn_comp(0.5)._ref_v, 'vpre', syn)
        ggn_vm_vec = h.Vector(ggn_vm[:, 1])
        tvec = h.Vector(ggn_vm[:, 0])
        model_dict[tvec] = None
        # vec.play(var_reference, t, continuous) for interpolating        
        ret = ggn_vm_vec.play(ggn_comp(0.5)._ref_v, tvec, 1)
        print('####', ret)
        model_dict[ggn_vm_vec] = None
        ggn_g_vec = h.Vector()
        ggn_g_vec.record(syn._ref_g)
        model_dict[ggn_g_vec] = None
    kc_vm_vec = h.Vector()
    kc_vm_vec.record(kc.soma(0.5)._ref_v)
    model_dict[kc_vm_vec] = None
    print('Built model')
    return (kc_vm_vec, ggn_g_vec)

def make_parser():
    parser = argparse.ArgumentParser(description='Simulate KC with GGN inhibition at multiple current injections')
    parser.add_argument('--kc-file', type=str, dest='kc_file',
                        required=True, help='KC cell template file'
                        ' (.hoc)')
    parser.add_argument('--kc', type=str, dest='kc', required=True,
                        help='KC cell template name in template file')
    parser.add_argument('--ggn-vm-file', type=str, action='append', dest='ggn_vm_file',
                        required=True,
                        help='CSV file with column 0 time in ms, column 1 GGN Vm in mV')
    parser.add_argument('--istart', type=str, dest='istart', help='Starting amplitude of current (with unit)')
    parser.add_argument('--iend', type=str, dest='iend', help='Ending amplitude of current (with unit)')
    parser.add_argument('--di', type=str, dest='di', help='Current increments (with unit)')
    parser.add_argument('--tstart', type=str, dest='tstart', help='Current injection start time (with unit)')
    parser.add_argument('--tend', type=str, dest='tend', help='Current injection end time (with unit)')    
    return parser


def main():
    parser = make_parser()
    args = parser.parse_args()
    logger.info('Command line args: {}'.format(str(sys.argv)))
    print(args.ggn_vm_file)
    # KCs with GGN inhibition
    inhibited_vec = defaultdict(list)
    solo_vec_list = []
    tstart = Q_(args.tstart).to('ms').m
    tend = Q_(args.tend).to('ms').m
    istart = Q_(args.istart).to('nA').m
    iend = Q_(args.iend).to('nA').m
    di = Q_(args.di).to('nA').m
    irange = np.arange(istart, iend + di/2.0, di)
    logger.info('Starting current: {} nA'.format(istart))
    logger.info('End current: {} nA'.format(iend))
    logger.info('Increment: {} nA'.format(di))
    logger.info('current range: {}'.format(irange))
    ggn_vm = {}
    for input_file in args.ggn_vm_file:
        ggn_vm[input_file] = np.loadtxt(input_file)
    for inject in irange:
        for input_file, vm in ggn_vm.items():
            kc_vvec, ggn_gvec = make_kc_with_dynaclamp(args.kc, args.kc_file, inject, tstart, tend, vm)
            inhibited_vec[input_file].append((kc_vvec, ggn_gvec))
        # KC without any inhibition
        kc_vvec, ggn_gvec = make_kc_with_dynaclamp(args.kc, args.kc_file, inject, tstart, tend)
        solo_vec_list.append(kc_vvec)
    tvec = h.Vector()
    tvec.record(h._ref_t)
    h.tstop = tend
    print('Init')    
    h.init()
    print('Run')
    h.run()
    print('Finished simulation')
    fig, ax = plt.subplots(nrows=len(irange)+1, ncols=len(ggn_vm)+1, sharex='all', sharey='all')
    t = np.array(tvec.x)
    solo_data = []
    for ii, vvec in enumerate(solo_vec_list):
        ax[ii+1, 0].plot(tvec, vvec, color='#e66101')
        solo_data.append(np.array(vvec.x))
    combined = np.vstack(solo_data)
    
    prefix = 'UTC' + timestamp.strftime('%Y%m%d_%H%M%S')
    fname = '{}_solo_kc.npz'.format(prefix)
    np.savez(fname,
             t=t,
             vm=combined,
             inject=irange)
    logger.info('Saved solo KC data in {}'.format(fname))
    for jj, input_file in enumerate(args.ggn_vm_file):
        fname = '{}_{}.npz'.format(prefix, os.path.basename(input_file))
        data = []
        kc_vm_list = inhibited_vec[input_file]
        for ii, (vvec, gvec) in enumerate(kc_vm_list):
            data.append(np.array(vvec.x))
            ax[ii+1, jj+1].plot(tvec, vvec, color='#e66101')
            ax[ii+1, 0].set_ylabel('{} pA'.format(irange[ii]*1e3))
            # ax[ii+1, 0].set_ylabel('{} pA'.format(int(np.round(irange[ii]*1e3))))  # to avoid decimal point when integer values 
        # ax[0, jj+1].plot(tvec, gvec)
        # ax[0, jj+1].plot(ggn_vm[input_file][:,0], ggn_vm[input_file][:,1])
        ax[0, jj+1].set_title(input_file)
        combined = np.vstack(data)
        np.savez(fname, combined=combined, irange=irange, ggn_vm=ggn_vm[input_file])
        logger.info('Saved data from dynamic clamp with input from {} in {}'.format(
            input_file, fname))
    for axis in ax.flat:
        axis.set_xlim(250, 1750)
        fig.set_size_inches(210/25.4, 290/25.4)
    fig.tight_layout()
    fig.savefig('{}_KC_dynamic_range_with_ggn_vm.svg'.format(prefix))
    plt.show()
    print('End')

    
if __name__ == '__main__':
    main()
    

# 
# kc_ggn_feedback_dclamp.py ends here
# kc_ggn_feedback_dclamp.py --- # Author: <NAME> # Created: Tue Aug 20 10:58:08 2019 (-0400) # Last-Updated: Wed Dec 11 17:32:49 2019 (-0500) #           By: <NAME> # Version: $Id$ # Code: This script for testing expansion of the dynamic range of a KC due to GGN inhibition.  Instead of running a whole simulation in the full network, we play the GGN membrane potential back to the KC.  We use GGN Vm from two simulations, one with low PN activity and another with high PN activity. # keep global reference of created model components so that they are # not garbage collected when out of scope Read KC model from `kc_file`, inject current `inject` nA, apply     dynamic clamp `ggn_vm`, which should be a 2D array with time (ms)     in column 0, and voltage (mV) in column 1. # vec.play(var_reference, t, continuous) for interpolating ###', ret) # KCs with GGN inhibition # KC without any inhibition # ax[ii+1, 0].set_ylabel('{} pA'.format(int(np.round(irange[ii]*1e3))))  # to avoid decimal point when integer values # ax[0, jj+1].plot(tvec, gvec) # ax[0, jj+1].plot(ggn_vm[input_file][:,0], ggn_vm[input_file][:,1]) # # kc_ggn_feedback_dclamp.py ends here