<jupyter_start><jupyter_text>Antennal neural network with intrinsic set-point (simulation)Here, we will simulate the antennal neural network based on both the anatomy and the control theoretic models. We will begin by simulating the network with an intrinsic set-point that arises due to the network connectivity. Later on (in another notebook) we will simulate it with varying levels of set-point inputs.<jupyter_code># Import necessary libraries
import nest
import numpy as np
from copy import deepcopy
from pprint import pprint
from collections import defaultdict
from pprint import pprint
#%matplotlib inline
%pylab inline<jupyter_output>Populating the interactive namespace from numpy and matplotlib<jupyter_text>Bristle field neural networkA class that defines the underlying connectivity of a single bristle field. Sensory neurons underlying the bristle field are either active or inactive, defined by the cuticle position. These neurons feed onto the antennal motor neurons (along with interneurons carrying set-point), which are simply integrate and fire neurons (iaf_psc_alpha). The motor neurons control muscle firing, which in turn is a low pass (exponential moving average) filter with a time constant of mTau. The muscle firing in turn indirectly controls cuticle position, based on other antennal muscle activity.<jupyter_code>class bbnn():
    '''The bbnn class defines the Bohm's bristle neural network. 
    It initializes sensory and motor neurons, and changes its rate 
    based on bristle activation. It saves all the data in a dict'''

    def __init__(self, sim_t, Bnum, delay, r_act, 
                 r_inact, r_in, weight, in_weight, 
                 mTau, CutPos):
        '''
        sim_t - Simulation step (default = 1 ms); 
        Bnum - Bristle number per field (default = 100 bristles); 
        delay - Synaptic delay in transmission (default = 1 ms);
        r_act - Activated sensory neuron firing rate (default = 50.0);
        r_inact - Inactivated sensory neuron firing rate (default = 10.0);
        r_in - Interneuron firing rate;
        weight - synaptic weight between sensory and motor neurons;
        in_weight - synaptic weight between interneuron and motor neurons;
        mTau - integration time of the muscle;
        CutPos - starting cuticle (/antenna) position;
        '''
        # Save constants
        self.sim_t = sim_t
        self.Bnum = Bnum
        self.weight = weight
        self.delay = delay
        self.in_weight = in_weight
        self.mTau = mTau
        
        # Initialize network
        self.initNetwork(CutPos, r_act, r_inact, r_in)
        

    def initNetwork(self, CutPos, r_act, r_inact, r_in):
        '''Set up the network in NEST
        CutPos - starting cuticle (/antenna) position;
        r_act - firing rate of activated bristles;
        r_inact - firing rate of inactivated bristles
        r_in - firing rate of set-point interneuron'''
        # Save act, inact rates
        self.r_act = r_act
        self.r_inact = r_inact
        self.r_in = r_in

        # Obtain bristle activity rate
        sn_rate = self.getFieldRate(CutPos, r_act, r_inact)

        # Declare Nodes
        ## Initialize poisson generators
        self.sn_pgs = nest.Create("poisson_generator", self.Bnum, sn_rate)
        self.jo = nest.Create("poisson_generator", 1,
                              {"rate": float(self.r_in)})
        ## Initialize neurons
        self.sn = nest.Create("parrot_neuron", self.Bnum)  # sensory neuron
        self.mn = nest.Create("iaf_psc_alpha")  # motor neuron
        self.ms_in = nest.Create("parrot_neuron", 1)  # interneuron
        ## Intialize spike detector
        self.sn_sd = nest.Create("spike_detector", 1)
        self.in_sd = nest.Create("spike_detector", 1)
        self.mn_sd = nest.Create("spike_detector", 1)

        # Setup Connections
        ## Poisson generators to parrot neurons
        nest.Connect(self.sn_pgs, self.sn, 'one_to_one')
        nest.Connect(self.jo, self.ms_in, 'one_to_one')
        ## Sensory to motor neurons
        nest.Connect(self.sn, self.mn, 'all_to_all', {
            "model": "static_synapse",
            "weight": self.weight,
            "delay": self.delay
        })
        nest.Connect(self.ms_in, self.mn, 'all_to_all', {
            "model": "static_synapse",
            "weight": self.in_weight * self.weight,
            "delay": self.delay
        })
        ## Neurons to spike detectors
        nest.Connect(self.sn, self.sn_sd)
        nest.Connect(self.ms_in, self.in_sd)
        nest.Connect(self.mn, self.mn_sd)

        nest.PrintNetwork()  # Not working for some reason
        # NEST network setup complete

        # Muscle activity related datastructure
        self.time = 0.0
        self.curr_ind = 0
        self.muscle_activity = []
        self.muscle_activity.append((self.time, 0.0))

        # Initialization complete

    def getMuscleActivity(self):
        ''' Obtain muscle firing rate based on the simulation time '''
        self.time += self.sim_t
        self.curr_ind += 1
        muscle_rate = moving_avg_rate(
            self.time,
            nest.GetStatus(self.mn_sd)[0]['events']['times'], self.mTau)
        self.muscle_activity.append((self.time, muscle_rate))
        return muscle_rate

    def changeRates(self, CutPos):
        ''' Change the firing rate of the poisson generators 
        based on bristle activity (cuticle position)'''
        # Get Rate
        rate = self.getFieldRate(CutPos, self.r_act, self.r_inact)
        # Set rate
        nest.SetStatus(self.sn_pgs, rate)

    def getFieldRate(self, CutPos, r_act, r_inact):
        ''' Obtain activity rate of each of the sensory neurons underneath the bristles '''
        Bfield_closed = np.concatenate(
            (np.zeros(self.Bnum - CutPos),
             np.ones(CutPos)))  # 1 to CutPos are active bristles
        Bfield_open = np.concatenate(
            (np.ones(self.Bnum - CutPos),
             np.zeros(CutPos)))  # CutPos to end are inactive bristles
        # Arbitary defintion of active and inactive. Simply change r_act, r_inact to inverse it
        field = ((Bfield_open * r_inact) + (Bfield_closed * r_act))

        return [{"rate": m} for m in field]

def moving_avg_rate(t, spike_times, tau):
    ''' Find rate at time t based on tau 
    assuming an moving average of spikes 
    with a window size of tau'''
     
    ind = len(spike_times) - 1
    muscle_rate = 0
    # Sum up spikes based on calcium decay
    while ((ind >= 0) and (t - spike_times[ind] < tau)):
        muscle_rate += 1000/tau
        ind -= 1
    
    return muscle_rate<jupyter_output><empty_output><jupyter_text>Antennal neural networkA overarching class representing the antennal neural network. This class creates and connects bristle field neural network objects, two in this simulation, and updates cuticle/antennal position based on muscle activity of each of the bristle field neural network.<jupyter_code>class ann():
    '''Antennal neural network class.'''

    def __init__(self, CutLoc, sim_t=10.0, Bnum=100, weight=50.0, 
                 delay=1.0, in_weight=1.0, mTau=50.0, mKp=50.0,
                 r_act=50.0, r_inact=5.0, r_in=50.0):
        '''Set defaults for the antennal neural network
        and initialize the network.
        
        Inputs:
            CutLoc - starting cuticle location'''

        # Set defaults
        self.sim_t = sim_t
        self.Bnum = Bnum
        self.weight = weight
        self.delay = delay
        self.in_weight = in_weight
        # muscle properties
        self.mTau = mTau
        self.mKp = mKp
        # neuron characteristics
        self.r_act = r_act
        self.r_inact = r_inact
        self.r_in = r_in
        # cuticle position
        self.CutLoc = CutLoc
        self.cuticle_position = []

        # Initialize field 1
        self.field1 = bbnn(self.sim_t, self.Bnum, self.delay, self.r_act,
                           self.r_inact, self.r_in, self.weight,
                           self.in_weight, self.mTau, self.CutLoc)

        # Initialize field 2 (flipped activation and inactivation rates
        self.field2 = bbnn(self.sim_t, self.Bnum, self.delay, self.r_inact,
                           self.r_act, self.r_in, self.weight, self.in_weight,
                           self.mTau, self.CutLoc)

    def Simulate(self, simulation_time):
        '''Simulate the neural network.'''
        time = np.arange(0, simulation_time, self.sim_t)
        CutLoc = self.CutLoc        

        # Simulate till end of time (literally :P)
        for t in time:
            # Simulate all networks
            nest.Simulate(self.sim_t)

            # Obtain muscle activity
            firing1 = self.field1.getMuscleActivity()
            firing2 = self.field2.getMuscleActivity()

            # Update position based on activity
            CutShift = np.round((firing1 - firing2) / self.mKp)
            CutLoc = int(CutLoc - CutShift)
            if CutLoc < 0: CutLoc = 0
            if CutLoc > self.Bnum: CutLoc = self.Bnum

            # Change rates
            self.field1.changeRates(CutLoc)
            self.field2.changeRates(CutLoc)

            # Append cuticle location
            self.cuticle_position.append((t, CutLoc))

            # Continue simulation
        # Simulation done; Extract data
        field1_activity = self.getBristleFieldActivity(self.field1)
        field2_activity = self.getBristleFieldActivity(self.field2)

        return field1_activity, field2_activity

    def getBristleFieldActivity(self, field):
        '''Extract relevant activity details from the field'''        
        
        # Extract sensory neuron spike trains
        sn_spikes = defaultdict(list)
        sn_spikes_raw = zip(
            nest.GetStatus(field.sn_sd)[0]['events']['times'],
            nest.GetStatus(field.sn_sd)[0]['events']['senders'])
        for time, sender in sn_spikes_raw:
            sn_spikes[sender].append(time)
        
        # Extract motor neuron spike trains
        mn_spikes = defaultdict(list)
        mn_spikes_raw = zip(
            nest.GetStatus(field.mn_sd)[0]['events']['times'],
            nest.GetStatus(field.mn_sd)[0]['events']['senders'])
        for time, sender in mn_spikes_raw:
            mn_spikes[sender].append(time)
        
        # Extract interneuron spike trains        
        in_spikes = defaultdict(list)
        in_spikes_raw = zip(
            nest.GetStatus(field.in_sd)[0]['events']['times'],
            nest.GetStatus(field.in_sd)[0]['events']['senders'])
        for time, sender in in_spikes_raw:
            in_spikes[sender].append(time)
            
        # Extract muscle activity
        muscle_activity = field.muscle_activity
        
        return (sn_spikes, mn_spikes, in_spikes, muscle_activity)
nest.ResetKernel()
network_1 = ann(0, sim_t=10.0, weight=50, mKp=25)
field1_activity, field2_activity = network_1.Simulate(500)

sn_spikes1, mn_spikes1, in_spikes1, muscle_activity1 = field1_activity 
muscle_activity1 = zip(*muscle_activity1)

sn_spikes2, mn_spikes2, in_spikes2, muscle_activity2 = field2_activity 
muscle_activity2 = zip(*muscle_activity2)

cuticle_position = zip(*network_1.cuticle_position)
subplot(311)
plot(cuticle_position[0], cuticle_position[1])
subplot(312)
plot(muscle_activity1[0], muscle_activity1[1])
plot(muscle_activity2[0], muscle_activity2[1])
subplot(313)
plot(muscle_activity2[0], np.array(muscle_activity2[1]) - np.array(muscle_activity1[1]))
nest.ResetKernel()
network_1 = ann(0, sim_t=1.0, mKp=125)
field1_activity, field2_activity = network_1.Simulate(500)

sn_spikes1, mn_spikes1, in_spikes1, muscle_activity1 = field1_activity 
muscle_activity1 = zip(*muscle_activity1)

sn_spikes2, mn_spikes2, in_spikes2, muscle_activity2 = field2_activity 
muscle_activity2 = zip(*muscle_activity2)

cuticle_position = zip(*network_1.cuticle_position)
subplot(311)
plot(cuticle_position[0], cuticle_position[1])
subplot(312)
plot(muscle_activity1[0], muscle_activity1[1])
plot(muscle_activity2[0], muscle_activity2[1])
subplot(313)
plot(muscle_activity2[0], np.array(muscle_activity2[1]) - np.array(muscle_activity1[1]))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09020344799411717, 'fraction_numerical': 0.015850968216357546, 'mean_word_length': 2.624518803671898, 'pattern_counts': {'":': 8, '<': 14, '<?xml version=': 0, '>': 14, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '491068', 'license': 'MIT', 'n_tokens_mistral': 3869, 'n_tokens_neox': 3692, 'n_words': 1001, 'path': 'Neural Circuit Simulation/.ipynb_checkpoints/nn-simulation-intrinsic-setpoint-rate-model-checkpoint.ipynb', 'repo_name': 'AbstractGeek/2019-Tuneable_reflexes_control_antennal_positioning_in_flying_hawkmoths', 'size': '12239'}