"""

Description: File listing procedures and functions.

Edit History: Created by Nilapratim Sengupta in July-August 2023.
              Edited by Nilapratim Sengupta in September-October 2023 for implementing parallelization.
              Edited by Nilapratim Sengupta in December 2023 to update options for location_start.
              Edited by Nilapratim Sengupta in January 2024 to make it generic with respect to axon/collateral diameter.
              Edited by Nilapratim Sengupta in January 2024 to update nseg values for the nodes and paranodes.
              Edited by Nilapratim Sengupta in May 2024 to calculate conduction time considering corresponding spikes.

"""

# Import statements
from neuron import h
from neuron.units import ms, mV, µm
from datetime import datetime
import math
import numpy as np
# import matplotlib.pyplot as plt
import pandas as pd

# Importing necessary files
import fileParameters
import fileCellCreation
import fileProceduresFunctions

# Global parameters
numOfCompartments_paranode = 8
numOfCompartments_juxtaparanode = 2
numOfCompartments_internode = 5

delta_paranode = 0.002
delta_internode = 0.02

specificCapacitance_axon = 1  # (uF/cm^2)
specificCapacitance_myelin = 0.5  # (uF/cm^2)

resistivity_axoplasm = 150
resistivity_periaxon = 150
resistivity_myelin = 300
resistivity_tightJunction = 600 # 600 Ohm.cm^2

# Initializing MPI and creating parallel context
h.nrnmpi_init()  
parallelContext = h.ParallelContext()





"""
Name: func_readParametersFromFile()
Type: Function

Input Parameter(s): 1. Filename
Return Parameter: Matrix of axon parameters 

Description: Function to read axon parameters for the population from a file.

Edit History: Created by Nilapratim Sengupta in August 2023.
"""
def func_readParametersFromFile(filename):
    
    # Initializing variables to track population size and parameter count
    size_population = 0
    count_axonParameters = 0
    list_tempLines = [] # List to temporarily store lines from the file

    # Opening the file for reading
    with open(filename, 'r') as file_axonParameters:
        for line in file_axonParameters:
            size_population += 1
            count_axonParameters = max(count_axonParameters, len(line.split('\t')))
            list_tempLines.append(line.strip()) # Storing stripped lines in the list

    # Creating numpy matrix to hold the axon parameters
    matrix_axonParameters = np.zeros((size_population, count_axonParameters))

    # Populating the matrix with axon parameters from the list
    for index, line in enumerate(list_tempLines):
        list_axonParameters = list(map(float, line.split('\t')))
        matrix_axonParameters[index, :] = list_axonParameters

    # Matrix to be returned
    return matrix_axonParameters

# End of func_readParametersFromFile()





"""
Name: func_readPerturbationsFromFile()
Type: Function

Input Parameter(s): 1. Filename
Return Parameter: Matrix of perturbations 

Description: Function to read perturbation schemes from a file.

Edit History: Created by Nilapratim Sengupta in August 2023.
"""
def func_readPerturbationsFromFile(filename):
    
    # Initializing variables to track number of segments and count of perturbation columns
    count_segments = 0
    count_perturbationColumns = 0
    list_tempLines = [] # List to temporarily store lines from the file

    # Opening the file for reading
    with open(filename, 'r') as file_perturbations:
        for line in file_perturbations:
            count_segments += 1
            count_perturbationColumns = max(count_perturbationColumns, len(line.split('\t')))
            list_tempLines.append(line.strip()) # Storing stripped lines in the list

    # Creating numpy matrix to hold the perturbation schemes
    matrix_perturbations = np.zeros((count_segments, count_perturbationColumns))

    # Populating the matrix with perturbation schemes from the list
    for index, line in enumerate(list_tempLines):
        list_perturbations = list(map(float, line.split('\t')))
        matrix_perturbations[index, :] = list_perturbations

    # Matrix to be returned
    return matrix_perturbations

# End of func_readPerturbationsFromFile()





"""
Name: func_updateConnectionSiteList_pyramidalCollateral()
Type: Function

Input Parameter(s): 1. Number of collaterals
                    2. List of connection site(s) along the axon for the collateral(s) 
                    3. Initial list of perturbation codes for the axon
Return Parameter: List (updated) of connection site(s) along the axon for the collateral(s) 

Description: Function to update the list of connection site(s) along the axon for the collateral(s).

Edit History: Created by Nilapratim Sengupta in September 2023.
"""
def func_updateConnectionSiteList_pyramidalCollateral(count_collaterals, connectionSiteList_collaterals, initialList_perturbationCode_axon):

    # Looping over the collateral(s)
    connectionSiteList_collaterals_updated = []
    for loopCounter_collateral in range(count_collaterals):
        connectionSite = connectionSiteList_collaterals[loopCounter_collateral]
        if (connectionSite > 0):
            count_extraNode = 0
            for loopCounter in range(len(initialList_perturbationCode_axon)):
                if (initialList_perturbationCode_axon[loopCounter] > 1 and loopCounter < connectionSite):
                    count_extraNode += (initialList_perturbationCode_axon[loopCounter] - 1) # If n shorter segments, n-1 new nodes
        connectionSiteList_collaterals_updated.append(int(connectionSite + count_extraNode))
        
    # Returning the updated list
    return connectionSiteList_collaterals_updated

# End of func_updateConnectionSiteList_pyramidalCollateral()





"""
Name: func_calculateSectionLengths()
Type: Function

Input Parameter(s): 1. Diameter of the process (axon/collateral) (um)
                    2. Length of each node um)
                    3. Lenth of each myelinated segment (um)
                    4. Perturbation code (0: unperturbed, 1: demyelinated, >=2: remyelinated)
                    5. Paranode-scaling code (0: fixed, 1: scaled)
Return Parameter: List of section lengths - [compartmentalLength_paranode, compartmentalLength_juxtaparanode, totalLength_internode]

Description: Function to calculate lengths of sections within myelinated segments.

Edit History: Created by Nilapratim Sengupta in August 2023.
              Edited by Nilapratim Sengupta in January 2024 to make it generic with respect to axon/collateral diameter.
"""
def func_calculateSectionLengths(diameter_process, length_node, length_segment, code_perturbation, code_paranodeScaling):


    # Lengths - compartmental and total
    if (code_perturbation < 2):
        compartmentalLength_paranode = (3.94 * (diameter_process / 4))
        totalLength_paranode = (numOfCompartments_paranode * compartmentalLength_paranode)
        length_extraNodes = 0
        compartmentalLength_juxtaparanode = (((length_segment - length_extraNodes) - totalLength_paranode) / (numOfCompartments_juxtaparanode + numOfCompartments_internode))
        totalLength_internode = (numOfCompartments_internode * compartmentalLength_juxtaparanode)
    elif (code_perturbation >= 2 and code_paranodeScaling == 0):
        compartmentalLength_paranode = (3.94 * (diameter_process / 4))
        totalLength_paranode = (numOfCompartments_paranode * compartmentalLength_paranode)
        length_extraNodes = (length_node * (code_perturbation - 1))
        compartmentalLength_juxtaparanode = (((length_segment - length_extraNodes) - (code_perturbation * totalLength_paranode)) / (code_perturbation * (numOfCompartments_juxtaparanode + numOfCompartments_internode)))
        totalLength_internode = (numOfCompartments_internode * compartmentalLength_juxtaparanode)
    elif (code_perturbation >= 2 and code_paranodeScaling == 1):
        compartmentalLength_paranode = ((3.94 * (diameter_process / 4)) / code_perturbation)
        totalLength_paranode = (numOfCompartments_paranode * compartmentalLength_paranode)
        length_extraNodes = (length_node * (code_perturbation - 1))
        compartmentalLength_juxtaparanode = (((length_segment - length_extraNodes) - (code_perturbation * totalLength_paranode)) / (code_perturbation * (numOfCompartments_juxtaparanode + numOfCompartments_internode)))
        totalLength_internode = (numOfCompartments_internode * compartmentalLength_juxtaparanode)
    
    # List to be returned
    list_sectionLengths = []
    list_sectionLengths.append(compartmentalLength_paranode)
    list_sectionLengths.append(compartmentalLength_juxtaparanode)
    list_sectionLengths.append(totalLength_internode)
    
    return list_sectionLengths

# End of func_calculateSectionLengths()
    
    
    


"""
Name: func_calculateExtracellularParameters()
Type: Function

Input Parameter(s): 1. Diameter of the process (axon/collateral) (um)
                    2. Length of each node um)
                    3. Lenth of each myelinated segment (um)
                    4. Number of myelin lamellae over each segment
                    5. Thickness of each myelin lamella (um)
                    6. Perturbation code (0: unperturbed, 1: demyelinated, >=2: remyelinated)
                    7. Paranode-scaling code (0: fixed, 1: scaled)
Return Parameter: List of lists (of parameter values for extracellular implementation) - [[xraxial0],[xraxial1],[xg0],[xg1],[xc0],[xc1],[e]]

Description: Function to calculate values of parameters for control/unperturbed myelinated segments.
             This function is based on calculations listed in an excel sheet, prepared by Gow-Devaux.

Edit History: Created by Nilapratim Sengupta in July 2023.
              Edited by Nilapratim Sengupta in January 2024 to make it generic with respect to axon/collateral diameter.
"""
def func_calculateExtracellularParameters(diameter_process, length_node, length_segment, count_lamellae, thickness_lamella, code_perturbation, code_paranodeScaling):

    # Lengths - compartmental and total
    [compartmentalLength_paranode, compartmentalLength_juxtaparanode, totalLength_internode] = func_calculateSectionLengths(diameter_process, length_node, length_segment, code_perturbation, code_paranodeScaling)
        
    # Cross-sectional areas
    crosssectionArea_node = (math.pi * (diameter_process / 2) ** 2) 
    crosssectionArea_paranodeOne = crosssectionArea_node
    crosssectionArea_paranodeTwo = crosssectionArea_node
    crosssectionArea_paranodeThree = crosssectionArea_node
    crosssectionArea_paranodeFour = crosssectionArea_node
    crosssectionArea_juxtaparanode = crosssectionArea_node
    crosssectionArea_internode = crosssectionArea_node

    # Peri-axonal annular areas
    periaxonAnnularArea_node = ((math.pi * ((diameter_process / 2) + delta_paranode) ** 2) - crosssectionArea_node)  
    periaxonAnnularArea_paranodeOne = ((math.pi * ((diameter_process / 2) + delta_paranode) ** 2) - crosssectionArea_paranodeOne) 
    periaxonAnnularArea_paranodeTwo = ((math.pi * ((diameter_process / 2) + delta_paranode) ** 2) - crosssectionArea_paranodeTwo) 
    periaxonAnnularArea_paranodeThree = ((math.pi * ((diameter_process / 2) + delta_paranode) ** 2) - crosssectionArea_paranodeThree) 
    periaxonAnnularArea_paranodeFour = ((math.pi * ((diameter_process / 2) + delta_paranode) ** 2) - crosssectionArea_paranodeFour) 
    periaxonAnnularArea_juxtaparanode = ((math.pi * ((diameter_process / 2) + delta_internode) ** 2) - crosssectionArea_juxtaparanode)  
    periaxonAnnularArea_internode = ((math.pi * ((diameter_process / 2) + delta_internode) ** 2) - crosssectionArea_internode)

    # Surface areas
    axonSurfaceArea_node = (math.pi * diameter_process * length_node)  
    axonSurfaceArea_paranodeOne = (math.pi * diameter_process * compartmentalLength_paranode)  
    axonSurfaceArea_paranodeTwo = (math.pi * diameter_process * compartmentalLength_paranode)  
    axonSurfaceArea_paranodeThree = (math.pi * diameter_process * compartmentalLength_paranode)  
    axonSurfaceArea_paranodeFour = (math.pi * diameter_process * compartmentalLength_paranode)  
    axonSurfaceArea_juxtaparanode = (math.pi * diameter_process * compartmentalLength_juxtaparanode) 
    axonSurfaceArea_internode = (math.pi * diameter_process * totalLength_internode) 

    # Myelin thicknesses
    myelinThickness_internode = (count_lamellae * thickness_lamella)
    myelinThickness_juxtaparanode = myelinThickness_internode
    myelinThickness_paranodeFour = (0.8 * myelinThickness_internode)
    myelinThickness_paranodeThree = (0.6 * myelinThickness_internode)
    myelinThickness_paranodeTwo = (0.4 * myelinThickness_internode)
    myelinThickness_paranodeOne = (0.2 * myelinThickness_internode)

    # Myelin lamellae
    myelinLamellae_internode = count_lamellae
    myelinLamellae_juxtaparanode = myelinLamellae_internode
    myelinLamellae_paranodeFour = (0.8 * myelinLamellae_internode)
    myelinLamellae_paranodeThree = (0.6 * myelinLamellae_internode)
    myelinLamellae_paranodeTwo = (0.4 * myelinLamellae_internode)
    myelinLamellae_paranodeOne = (0.2 * myelinLamellae_internode)

    # Myelin surface areas
    myelinSurfaceArea_paranodeOne = math.pi * compartmentalLength_paranode * (diameter_process + (2 * delta_paranode) + (2 * myelinThickness_paranodeOne)) 
    myelinSurfaceArea_paranodeTwo = math.pi * compartmentalLength_paranode * (diameter_process + (2 * delta_paranode) + (2 * myelinThickness_paranodeTwo))  
    myelinSurfaceArea_paranodeThree = math.pi * compartmentalLength_paranode * (diameter_process + (2 * delta_paranode) + (2 * myelinThickness_paranodeThree)) 
    myelinSurfaceArea_paranodeFour = math.pi * compartmentalLength_paranode * (diameter_process + (2 * delta_paranode) + (2 * myelinThickness_paranodeFour))  
    myelinSurfaceArea_juxtaparanode = math.pi * compartmentalLength_juxtaparanode * (diameter_process + (2 * delta_internode) + (2 * myelinThickness_juxtaparanode))  
    myelinSurfaceArea_internode = math.pi * totalLength_internode * (diameter_process + (2 * delta_internode) + (2 * myelinThickness_internode))  

    # Axoplasmic resistances
    axoplasmResistance_node = resistivity_axoplasm * length_node / crosssectionArea_node  # (10^4 ohm)
    axoplasmResistance_paranodeOne = resistivity_axoplasm * compartmentalLength_paranode / crosssectionArea_paranodeOne  # (10^4 ohm)
    axoplasmResistance_paranodeTwo = resistivity_axoplasm * compartmentalLength_paranode / crosssectionArea_paranodeTwo  # (10^4 ohm)
    axoplasmResistance_paranodeThree = resistivity_axoplasm * compartmentalLength_paranode / crosssectionArea_paranodeThree  # (10^4 ohm)
    axoplasmResistance_paranodeFour = resistivity_axoplasm * compartmentalLength_paranode / crosssectionArea_paranodeFour  # (10^4 ohm)
    axoplasmResistance_juxtaparanode = resistivity_axoplasm * compartmentalLength_juxtaparanode / crosssectionArea_juxtaparanode  # (10^4 ohm)
    axoplasmResistance_internode = resistivity_axoplasm * totalLength_internode / crosssectionArea_internode  # (10^4 ohm)

    # Peri-axonal resistances
    periaxonResistance_paranodeOne = resistivity_periaxon * compartmentalLength_paranode / periaxonAnnularArea_paranodeOne  # (10^4 ohm)
    periaxonResistance_paranodeTwo = resistivity_periaxon * compartmentalLength_paranode / periaxonAnnularArea_paranodeTwo  # (10^4 ohm)
    periaxonResistance_paranodeThree = resistivity_periaxon * compartmentalLength_paranode / periaxonAnnularArea_paranodeThree  # (10^4 ohm)
    periaxonResistance_paranodeFour = resistivity_periaxon * compartmentalLength_paranode / periaxonAnnularArea_paranodeFour  # (10^4 ohm)
    periaxonResistance_juxtaparanode = resistivity_periaxon * compartmentalLength_juxtaparanode / periaxonAnnularArea_juxtaparanode  # (10^4 ohm)
    periaxonResistance_internode = resistivity_periaxon * totalLength_internode / periaxonAnnularArea_internode  # (10^4 ohm)

    # Axonal capacitances
    axonCapacitance_node = specificCapacitance_axon * axonSurfaceArea_node  # (10^-8 uF)
    axonCapacitance_paranodeOne = specificCapacitance_axon * axonSurfaceArea_paranodeOne  # (10^-8 uF)
    axonCapacitance_paranodeTwo = specificCapacitance_axon * axonSurfaceArea_paranodeTwo  # (10^-8 uF)
    axonCapacitance_paranodeThree = specificCapacitance_axon * axonSurfaceArea_paranodeThree  # (10^-8 uF)
    axonCapacitance_paranodeFour = specificCapacitance_axon * axonSurfaceArea_paranodeFour  # (10^-8 uF)
    axonCapacitance_juxtaparanode = specificCapacitance_axon * axonSurfaceArea_juxtaparanode  # (10^-8 uF)
    axonCapacitance_internode = specificCapacitance_axon * axonSurfaceArea_internode  # (10^-8 uF)

    # Myelin capacitances
    myelinCapacitance_paranodeOne = specificCapacitance_myelin * myelinSurfaceArea_paranodeOne / myelinLamellae_paranodeOne  # (10^-8 uF)
    myelinCapacitance_paranodeTwo = specificCapacitance_myelin * myelinSurfaceArea_paranodeTwo / myelinLamellae_paranodeTwo  # (10^-8 uF)
    myelinCapacitance_paranodeThree = (specificCapacitance_myelin * myelinSurfaceArea_paranodeThree / myelinLamellae_paranodeThree) # (10^-8 uF)
    myelinCapacitance_paranodeFour = (specificCapacitance_myelin * myelinSurfaceArea_paranodeFour / myelinLamellae_paranodeFour) # (10^-8 uF)
    myelinCapacitance_juxtaparanode = (specificCapacitance_myelin * myelinSurfaceArea_juxtaparanode / myelinLamellae_juxtaparanode) # (10^-8 uF)
    myelinCapacitance_internode = (specificCapacitance_myelin * myelinSurfaceArea_internode / myelinLamellae_internode) # (10^-8 uF)

    #  Myelin resistances
    myelinResistance_paranodeOne = (resistivity_myelin * myelinLamellae_paranodeOne / myelinSurfaceArea_paranodeOne) # (10^8 ohm) 
    myelinResistance_paranodeTwo = (resistivity_myelin * myelinLamellae_paranodeTwo / myelinSurfaceArea_paranodeTwo) # (10^8 ohm) 
    myelinResistance_paranodeThree = (resistivity_myelin * myelinLamellae_paranodeThree / myelinSurfaceArea_paranodeThree) # (10^8 ohm) 
    myelinResistance_paranodeFour = (resistivity_myelin * myelinLamellae_paranodeFour / myelinSurfaceArea_paranodeFour) # (10^8 ohm) 
    myelinResistance_juxtaparanode = (resistivity_myelin * myelinLamellae_juxtaparanode / myelinSurfaceArea_juxtaparanode) # (10^8 ohm) 
    myelinResistance_internode = (resistivity_myelin * myelinLamellae_internode / myelinSurfaceArea_internode) # (10^8 ohm)

    # Tight junctional resistances
    tightjunctionResistance_paranodeOne = (resistivity_tightJunction / (math.pi * (diameter_process + (2 * delta_paranode)) * compartmentalLength_paranode)) # (10^8 ohm)
    tightjunctionResistance_paranodeTwo = (resistivity_tightJunction / (math.pi * (diameter_process + (2 * delta_paranode)) * compartmentalLength_paranode)) # (10^8 ohm)
    tightjunctionResistance_paranodeThree = (resistivity_tightJunction / (math.pi * (diameter_process + (2 * delta_paranode)) * compartmentalLength_paranode)) # (10^8 ohm)
    tightjunctionResistance_paranodeFour = (resistivity_tightJunction / (math.pi * (diameter_process + (2 * delta_paranode)) * compartmentalLength_paranode)) # (10^8 ohm)
    tightjunctionResistance_juxtaparanode = (resistivity_tightJunction / (math.pi * (diameter_process + (2 * delta_internode)) * compartmentalLength_juxtaparanode)) # (10^8 ohm)
    tightjunctionResistance_internode = (resistivity_tightJunction / (math.pi * (diameter_process + (2 * delta_internode)) * totalLength_internode)) # (10^8 ohm)

    # xraxial[0] for NEURON
    xraxial0_node = (resistivity_periaxon / (periaxonAnnularArea_node * 0.00000001) * 0.000001) # (Mohm/cm)
    xraxial0_paranodeOne = (2 * (resistivity_periaxon / (periaxonAnnularArea_paranodeOne * 0.00000001) * 0.000001)) # (Mohm/cm)
    xraxial0_paranodeTwo = (2 * (resistivity_periaxon / (periaxonAnnularArea_paranodeTwo * 0.00000001) * 0.000001)) # (Mohm/cm)
    xraxial0_paranodeThree = (2 * (resistivity_periaxon / (periaxonAnnularArea_paranodeThree * 0.00000001) * 0.000001)) # (Mohm/cm)
    xraxial0_paranodeFour = (2 * (resistivity_periaxon / (periaxonAnnularArea_paranodeFour * 0.00000001) * 0.000001)) # (Mohm/cm)
    xraxial0_juxtaparanode = (resistivity_periaxon / (periaxonAnnularArea_juxtaparanode * 0.00000001) * 0.000001) # (Mohm/cm)
    xraxial0_internode = (resistivity_periaxon / (periaxonAnnularArea_internode * 0.00000001) * 0.000001) # (Mohm/cm)
    
    xraxial0_list = []
    xraxial0_list.append(xraxial0_node)
    xraxial0_list.append(xraxial0_paranodeOne)
    xraxial0_list.append(xraxial0_paranodeTwo)
    xraxial0_list.append(xraxial0_paranodeThree)
    xraxial0_list.append(xraxial0_paranodeFour)
    xraxial0_list.append(xraxial0_juxtaparanode)
    xraxial0_list.append(xraxial0_internode)

    # xraxial[1] for NEURON
    xraxial1_value = 1000000000
    xraxial1_node = xraxial1_value # (Mohm/cm)
    xraxial1_paranodeOne = xraxial1_value # (Mohm/cm)
    xraxial1_paranodeTwo = xraxial1_value # (Mohm/cm)
    xraxial1_paranodeThree = xraxial1_value # (Mohm/cm)
    xraxial1_paranodeFour = xraxial1_value # (Mohm/cm)
    xraxial1_juxtaparanode = xraxial1_value # (Mohm/cm)
    xraxial1_internode = xraxial1_value # (Mohm/cm)
    
    xraxial1_list = []
    xraxial1_list.append(xraxial1_node)
    xraxial1_list.append(xraxial1_paranodeOne)
    xraxial1_list.append(xraxial1_paranodeTwo)
    xraxial1_list.append(xraxial1_paranodeThree)
    xraxial1_list.append(xraxial1_paranodeFour)
    xraxial1_list.append(xraxial1_juxtaparanode)
    xraxial1_list.append(xraxial1_internode)

    # xg[0] for NEURON
    xg0_node = 1000000000 # (S/cm^2)
    xg0_paranodeOne = (1 / (tightjunctionResistance_paranodeOne * axonSurfaceArea_paranodeOne)) # (S/cm^2)
    xg0_paranodeTwo = (1 / (tightjunctionResistance_paranodeTwo * axonSurfaceArea_paranodeTwo)) # (S/cm^2)
    xg0_paranodeThree = (1 / (tightjunctionResistance_paranodeThree * axonSurfaceArea_paranodeThree)) # (S/cm^2)
    xg0_paranodeFour = (1 / (tightjunctionResistance_paranodeFour * axonSurfaceArea_paranodeFour)) # (S/cm^2)
    xg0_juxtaparanode = (1 / (tightjunctionResistance_juxtaparanode * axonSurfaceArea_juxtaparanode)) # (S/cm^2)
    xg0_internode = (1 / (tightjunctionResistance_internode * axonSurfaceArea_internode)) # (S/cm^2)
    
    xg0_list = []
    xg0_list.append(xg0_node)
    xg0_list.append(xg0_paranodeOne)
    xg0_list.append(xg0_paranodeTwo)
    xg0_list.append(xg0_paranodeThree)
    xg0_list.append(xg0_paranodeFour)
    xg0_list.append(xg0_juxtaparanode)
    xg0_list.append(xg0_internode)

    # xg[1] for NEURON
    xg1_node = 1000000000 # (S/cm^2)
    xg1_paranodeOne = (1 / (myelinResistance_paranodeOne * axonSurfaceArea_paranodeOne)) # (S/cm^2)
    xg1_paranodeTwo = (1 / (myelinResistance_paranodeTwo * axonSurfaceArea_paranodeTwo)) # (S/cm^2)
    xg1_paranodeThree = (1 / (myelinResistance_paranodeThree * axonSurfaceArea_paranodeThree)) # (S/cm^2)
    xg1_paranodeFour = (1 / (myelinResistance_paranodeFour * axonSurfaceArea_paranodeFour)) # (S/cm^2)
    xg1_juxtaparanode = (1 / (myelinResistance_juxtaparanode * axonSurfaceArea_juxtaparanode)) # (S/cm^2)
    xg1_internode = (1 / (myelinResistance_internode * axonSurfaceArea_internode)) # (S/cm^2)
    
    xg1_list = []
    xg1_list.append(xg1_node)
    xg1_list.append(xg1_paranodeOne)
    xg1_list.append(xg1_paranodeTwo)
    xg1_list.append(xg1_paranodeThree)
    xg1_list.append(xg1_paranodeFour)
    xg1_list.append(xg1_juxtaparanode)
    xg1_list.append(xg1_internode)

    # xc[0] for NEURON
    xc0_value = 0
    xc0_node = xc0_value # (uF/cm^2)
    xc0_paranodeOne = xc0_value # (uF/cm^2)
    xc0_paranodeTwo = xc0_value # (uF/cm^2)
    xc0_paranodeThree = xc0_value # (uF/cm^2)
    xc0_paranodeFour = xc0_value # (uF/cm^2)
    xc0_juxtaparanode = xc0_value # (uF/cm^2)
    xc0_internode = xc0_value # (uF/cm^2)
    
    xc0_list = []
    xc0_list.append(xc0_node)
    xc0_list.append(xc0_paranodeOne)
    xc0_list.append(xc0_paranodeTwo)
    xc0_list.append(xc0_paranodeThree)
    xc0_list.append(xc0_paranodeFour)
    xc0_list.append(xc0_juxtaparanode)
    xc0_list.append(xc0_internode)

    # xc[1] for NEURON
    xc1_node = 0 # (uF/cm^2)
    xc1_paranodeOne = (myelinCapacitance_paranodeOne / axonSurfaceArea_paranodeOne) # (uF/cm^2)
    xc1_paranodeTwo = (myelinCapacitance_paranodeTwo / axonSurfaceArea_paranodeTwo) # (uF/cm^2)
    xc1_paranodeThree = (myelinCapacitance_paranodeThree / axonSurfaceArea_paranodeThree) # (uF/cm^2)
    xc1_paranodeFour = (myelinCapacitance_paranodeFour / axonSurfaceArea_paranodeFour) # (uF/cm^2)
    xc1_juxtaparanode = (myelinCapacitance_juxtaparanode / axonSurfaceArea_juxtaparanode) # (uF/cm^2)
    xc1_internode = (myelinCapacitance_internode / axonSurfaceArea_internode) # (uF/cm^2)
    
    xc1_list = []
    xc1_list.append(xc1_node)
    xc1_list.append(xc1_paranodeOne)
    xc1_list.append(xc1_paranodeTwo)
    xc1_list.append(xc1_paranodeThree)
    xc1_list.append(xc1_paranodeFour)
    xc1_list.append(xc1_juxtaparanode)
    xc1_list.append(xc1_internode)

    # e_extracellular for NEURON
    eExtracellular_value = 0
    eExtracellular_node = eExtracellular_value # (mV)
    eExtracellular_paranodeOne = eExtracellular_value # (mV)
    eExtracellular_paranodeTwo = eExtracellular_value # (mV)
    eExtracellular_paranodeThree = eExtracellular_value # (mV)
    eExtracellular_paranodeFour = eExtracellular_value # (mV)
    eExtracellular_juxtaparanode = eExtracellular_value # (mV)
    eExtracellular_internode = eExtracellular_value # (mV)
    
    eExtracellular_list = []
    eExtracellular_list.append(eExtracellular_node)
    eExtracellular_list.append(eExtracellular_paranodeOne)
    eExtracellular_list.append(eExtracellular_paranodeTwo)
    eExtracellular_list.append(eExtracellular_paranodeThree)
    eExtracellular_list.append(eExtracellular_paranodeFour)
    eExtracellular_list.append(eExtracellular_juxtaparanode)
    eExtracellular_list.append(eExtracellular_internode)
    
    # List to be returned
    list_extracellularParameters = []
    list_extracellularParameters.append(xraxial0_list)
    list_extracellularParameters.append(xraxial1_list)
    list_extracellularParameters.append(xg0_list)
    list_extracellularParameters.append(xg1_list)
    list_extracellularParameters.append(xc0_list)
    list_extracellularParameters.append(xc1_list)
    list_extracellularParameters.append(eExtracellular_list)
    
    return list_extracellularParameters

# End of func_calculateExtracellularParameters()





"""
Name: proc_setGeometry_pyramidalSomaDendrites()
Type: Procedure

Input Parameter(s): 1. Object - pyramidalCell
                    2. Type of cell (1: default, 2: young; 3: aged)
Return Parameter: None

Description: Procedure to set length, diameter and nseg values for the soma and dendrites of the pyramidal cell.
             The corresponding values for the soma and the dendrites are from the Rumbell/Traub models.
                 
Edit History: Created by Nilapratim Sengupta in July-August 2023.
"""
def proc_setGeometry_pyramidalSomaDendrites(pyramidalCell, type_cell):
    
    # Defining geometry of the soma and dendrites
    for sec in pyramidalCell.pyramidalSomaDendrites.Soma:
        sec.L = 15
        sec.diam = 16
        
    for sec in pyramidalCell.pyramidalSomaDendrites.Dendrites:
        sec.L = 50
        
    for sec in pyramidalCell.pyramidalSomaDendrites.Distal:
        sec.diam = 1.6
        
    for sec in pyramidalCell.pyramidalSomaDendrites.Oblique:
        sec.diam = 1
        
    for sec in pyramidalCell.pyramidalSomaDendrites.Basal:
        sec.diam = 1
        
    for sec in pyramidalCell.pyramidalSomaDendrites.Aux:
        sec.L = 15 / 2
        sec.diam = 16
        
    for loopCounter in range(4):
        pyramidalCell.pyramidalSomaDendrites.aux10to13[loopCounter].L = 50 / 2
        pyramidalCell.pyramidalSomaDendrites.aux10to13[loopCounter].diam = 8
        
    pyramidalCell.pyramidalSomaDendrites.comp[38].diam = 8
    pyramidalCell.pyramidalSomaDendrites.comp[39].diam = 8 * 0.9
    pyramidalCell.pyramidalSomaDendrites.comp[40].diam = 8 * 0.8
    
    for sec in pyramidalCell.pyramidalSomaDendrites.Level8:
        sec.diam = 4
        
    for sec in pyramidalCell.pyramidalSomaDendrites.Level9:
        sec.diam = 4
        
    # Scaling of length and diameter of the soma-dendrite compartments based on cell type
    scaling_factors = {1: 0.888, 2: 0.724, 3: 0.794} # Dictionary of scaling_factors
    scaling_factor = scaling_factors.get(type_cell, 0.888) # Default scaling_factor = 0.888
    
    for sec in pyramidalCell.pyramidalSomaDendrites.SomaDendrites:
        sec.L *= scaling_factor
        sec.diam *= scaling_factor
            
# End of proc_setGeometry_pyramidalSomaDendrites()
            
            
 
            
            
"""
Name: proc_setGeometry_pyramidalAxon()
Type: Procedure

Input Parameter(s): 1. Object - pyramidalCell
                    2. Diameter of the axon (um)
                    3. Length of each node um)
                    4. Lenth of each myelinated segment (um)
                    5. Perturbation code (0: unperturbed, 1: demyelinated, >=2: remyelinated)
                    6. Paranode-scaling code (0: fixed, 1: scaled)
Return Parameter: None

Description: Procedure to set length, diameter and nseg values for different sections of the axon of the pyramidal cell.
             The corresponding values for the axon sections are calculated based on the Scurfield-Latimer/Gow-Devaux models.
                 
Edit History: Created by Nilapratim Sengupta in July-August 2023.
              Edited by Nilapratim Sengupta in January 2024 to update certain nseg values (as advised by Ted).
"""
def proc_setGeometry_pyramidalAxon(pyramidalCell, diameter_axon, length_node, length_segment, code_perturbation, code_paranodeScaling):
    
    # Defining geometry of the axon hillock and initial segment
    area_soma = 0
    for seg in pyramidalCell.pyramidalSomaDendrites.comp[1]:
        area_soma += seg.area()
        
    equivalentDiameter = (math.sqrt(area_soma / (4 * math.pi)) / 10)

    pyramidalCell.pyramidalAxon.hill_axon.L = 3 # 10 # Changed on 09/10/2024
    pyramidalCell.pyramidalAxon.hill_axon(0).diam = (4 * equivalentDiameter)
    pyramidalCell.pyramidalAxon.hill_axon(1).diam = equivalentDiameter
    pyramidalCell.pyramidalAxon.hill_axon.nseg = 5
    
    pyramidalCell.pyramidalAxon.iseg_axon.L = 30 # 15 # Changed on 09/10/2024
    pyramidalCell.pyramidalAxon.iseg_axon(0).diam = equivalentDiameter
    pyramidalCell.pyramidalAxon.iseg_axon(1).diam = diameter_axon
    pyramidalCell.pyramidalAxon.iseg_axon.nseg = 5
    
    # Obtaining lengths of the different compartments of the myelinated segments
    [compartmentalLength_paranode, compartmentalLength_juxtaparanode, totalLength_internode] = func_calculateSectionLengths(diameter_axon, length_node, length_segment, code_perturbation, code_paranodeScaling)
     
    # Defining geometry of the nodes and compartments of the myelinated segments along the axon
    for sec in pyramidalCell.pyramidalAxon.Nodes_axon:
        sec.L = length_node
        sec.diam = diameter_axon
        sec.nseg = 1 # 13
        
    for sec in pyramidalCell.pyramidalAxon.ParanodeOnes_axon:
        sec.L = compartmentalLength_paranode
        sec.diam = diameter_axon
        sec.nseg = 1 # 5
        
    for sec in pyramidalCell.pyramidalAxon.ParanodeTwos_axon:
        sec.L = compartmentalLength_paranode
        sec.diam = diameter_axon
        sec.nseg = 1 # 5
        
    for sec in pyramidalCell.pyramidalAxon.ParanodeThrees_axon:
        sec.L = compartmentalLength_paranode
        sec.diam = diameter_axon
        sec.nseg = 1 # 5
        
    for sec in pyramidalCell.pyramidalAxon.ParanodeFours_axon:
        sec.L = compartmentalLength_paranode
        sec.diam = diameter_axon
        sec.nseg = 1 # 5
        
    for sec in pyramidalCell.pyramidalAxon.Juxtaparanodes_axon:
        sec.L = compartmentalLength_juxtaparanode
        sec.diam = diameter_axon
        sec.nseg = 5
        
    for sec in pyramidalCell.pyramidalAxon.Internodes_axon:
        sec.L = totalLength_internode
        sec.diam = diameter_axon
        sec.nseg = 9
            
# End of proc_setGeometry_pyramidalAxon()





"""
Name: proc_setGeometry_pyramidalCollateral()
Type: Procedure

Input Parameter(s): 1. Object - pyramidalCell
                    2. Index of collateral
                    3. Diameter of the collateral (um)
                    4. Length of each node um)
                    5. Lenth of each myelinated segment (um)
                    6. Perturbation code (0: unperturbed, 1: demyelinated, >=2: remyelinated)
                    7. Paranode-scaling code (0: fixed, 1: scaled)
Return Parameter: None

Description: Procedure to set length, diameter and nseg values for different sections of the collateral(s) of the pyramidal cell.
             The corresponding values for the collateral sections are calculated based on the Scurfield-Latimer/Gow-Devaux models.
                 
Edit History: Created by Nilapratim Sengupta in July-August 2023.
              Edited by Nilapratim Sengupta in January 2024 to update collateral diameter.
              Edited by Nilapratim Sengupta in January 2024 to update certain nseg values (as advised by Ted).
              Edited by Nilapratim Sengupta in September 2024 to ensure collaterals can be modified individually.
"""
def proc_setGeometry_pyramidalCollateral(pyramidalCell, index_collateral, diameter_collateral, length_node, length_segment, code_perturbation, code_paranodeScaling):
    
    # Obtaining lengths of the different compartments of the myelinated segments
    [compartmentalLength_paranode, compartmentalLength_juxtaparanode, totalLength_internode] = func_calculateSectionLengths(diameter_collateral, length_node, length_segment, code_perturbation, code_paranodeScaling)
        
    for sec in pyramidalCell.pyramidalCollateral[index_collateral].Nodes_collateral:
        sec.L = length_node
        sec.diam = diameter_collateral
        sec.nseg = 1 # 13
        
    for sec in pyramidalCell.pyramidalCollateral[index_collateral].ParanodeOnes_collateral:
        sec.L = compartmentalLength_paranode
        sec.diam = diameter_collateral
        sec.nseg = 1 # 5
        
    for sec in pyramidalCell.pyramidalCollateral[index_collateral].ParanodeTwos_collateral:
        sec.L = compartmentalLength_paranode
        sec.diam = diameter_collateral
        sec.nseg = 1 # 5
        
    for sec in pyramidalCell.pyramidalCollateral[index_collateral].ParanodeThrees_collateral:
        sec.L = compartmentalLength_paranode
        sec.diam = diameter_collateral
        sec.nseg = 1 # 5
        
    for sec in pyramidalCell.pyramidalCollateral[index_collateral].ParanodeFours_collateral:
        sec.L = compartmentalLength_paranode
        sec.diam = diameter_collateral
        sec.nseg = 1 # 5
        
    for sec in pyramidalCell.pyramidalCollateral[index_collateral].Juxtaparanodes_collateral:
        sec.L = compartmentalLength_juxtaparanode
        sec.diam = diameter_collateral
        sec.nseg = 5
        
    for sec in pyramidalCell.pyramidalCollateral[index_collateral].Internodes_collateral:
        sec.L = totalLength_internode
        sec.diam = diameter_collateral
        sec.nseg = 9
            
# End of proc_setGeometry_pyramidalCollateral()





"""
Name: proc_setBiophysics_pyramidalSomaDendrites()
Type: Procedure

Input Parameter(s): 1. Object - pyramidalCell
                    2. dnap
                    3. dkc
Return Parameter: None

Description: Procedure to set active and passive membrane mechanism properties for the soma and dendrites of the pyramidal cell.
             The corresponding values for the soma and the dendrites are from the Rumbell/Traub models.
                 
Edit History: Created by Nilapratim Sengupta in July-August 2023.
"""
def proc_setBiophysics_pyramidalSomaDendrites(pyramidalCell):
    
    # Default values of specific parameters
    dnap = 0
    dkc = 1.6
    
    # Defining passive properties of the soma and dendrites
    for sec in pyramidalCell.pyramidalSomaDendrites.Soma:
        sec.Ra = 250
        sec.cm = 0.9
        sec.insert('pas')
        sec.g_pas = 2e-05
        sec.e_pas = -70
    
    for sec in pyramidalCell.pyramidalSomaDendrites.Dendrites:
        sec.Ra = 250
        sec.cm = 0.9
        sec.insert('pas')
        sec.g_pas = 2e-05
        sec.e_pas = -70
    
    for sec in pyramidalCell.pyramidalSomaDendrites.Aux:
        sec.Ra = 250
        sec.cm = 0
    
    # Inserting active ion channels in the soma and dendrites
    for sec in pyramidalCell.pyramidalSomaDendrites.Soma:
        sec.insert('cad')
        sec.insert('naf')
        sec.insert('nap')
        sec.insert('kdr')
        sec.insert('ka')
        sec.insert('kc')
        sec.insert('kahp') 
        sec.insert('k2')
        sec.insert('km')
        sec.insert('cat')
        sec.insert('cal')
        sec.insert('ar') 
        

    
    for sec in pyramidalCell.pyramidalSomaDendrites.Dendrites:
        sec.insert('cad')
        sec.insert('naf')
        sec.insert('nap')
        sec.insert('kdr')
        sec.insert('ka')
        sec.insert('kc')
        sec.insert('kahp') 
        sec.insert('k2')
        sec.insert('km')
        sec.insert('cat')
        sec.insert('cal')
        sec.insert('ar')
    
    # Setting ceiling value for calcium dynamics    
    h.ceiling_cad = 1000
    
    # Defining active properties of the soma and dendrites
    for sec in pyramidalCell.pyramidalSomaDendrites.Soma:
        sec.phi_cad = 52 / 2e-3
        sec.beta_cad = 1 / 100  # In the paper beta = 50 [ms]
        
        sec.gbar_naf = 150e-3 * 1.25
        sec.gbar_nap = dnap * 0.0032 * sec.gbar_naf
        sec.gbar_kdr = 125e-3
        sec.gbar_ka = 30e-3
        sec.gbar_kc = dkc * 7.5e-3  # In tha paper 'dkc * 12e-3'
        sec.gbar_kahp = 0.1e-3
        sec.gbar_k2 = 0.1e-3
        sec.gbar_km = 2.5 * 1.5e-3 * 2
        sec.gbar_cat = 0.1e-3
        sec.gbar_cal = 0.5e-3
        sec.gbar_ar = 0.25e-3

    for sec in pyramidalCell.pyramidalSomaDendrites.Dendrites:
        sec.phi_cad = 52 / 2e-3
        sec.beta_cad = 1 / 20
        
        sec.gbar_naf = 6.25e-3
        sec.gbar_nap = dnap * 0.0032 * sec.gbar_naf
        sec.gbar_kdr = 0
        sec.gbar_ka = 2e-3
        sec.gbar_kc = 0
        sec.gbar_kahp = 0.1e-3
        sec.gbar_k2 = 0.1e-3
        sec.gbar_km = 2.5 * 1.5e-3 * 2
        sec.gbar_cat = 0.1e-3
        sec.gbar_cal = 0.5e-3
        sec.gbar_ar = 0.25e-3

    for sec in pyramidalCell.pyramidalSomaDendrites.Proximal:
        sec.gbar_naf = 75e-3 * 1.25
        sec.gbar_nap = dnap * 0.0032 * sec.gbar_naf
        sec.gbar_kdr = 75e-3 * 1.25
        sec.gbar_kc = dkc * 7.5e-3  # In the paper 'dkc * 12e-3'

    for sec in pyramidalCell.pyramidalSomaDendrites.Distal:
        sec.gbar_cal = 3e-3

    pyramidalCell.pyramidalSomaDendrites.comp[38].gbar_ka = 30e-3
    pyramidalCell.pyramidalSomaDendrites.comp[38].gbar_naf = 125e-3
    pyramidalCell.pyramidalSomaDendrites.comp[38].gbar_nap = dnap * 0.0032 * pyramidalCell.pyramidalSomaDendrites.comp[38].gbar_naf  # In the FORTRAN code 0.004
    pyramidalCell.pyramidalSomaDendrites.comp[38].gbar_kdr = 125e-3  # In the paper '75e-3 * 1.25'
    pyramidalCell.pyramidalSomaDendrites.comp[38].gbar_kc = dkc * 7.5e-3  # In the paper 'dkc * 12e-3'

    for sec in pyramidalCell.pyramidalSomaDendrites.Soma:
        sec.ena = 50
        sec.ek = -95
        sec.eca = 125
        
    for sec in pyramidalCell.pyramidalSomaDendrites.Dendrites:
        sec.ena = 50
        sec.ek = -95
        sec.eca = 125

# End of proc_setBiophysics_pyramidalSomaDendrites()





"""
Name: proc_setBiophysics_pyramidalAxon()
Type: Procedure

Input Parameter(s): 1. Object - pyramidalCell
                    2. List of lists (of parameter values for extracellular implementation) - [[xraxial0],[xraxial1],[xg0],[xg1],[xc0],[xc1],[e]]
Return Parameter: None

Description: Procedure to set active and passive membrane mechanism properties for different sections of the axon of the pyramidal cell.
             The corresponding values for the axon sections are calculated based on the Scurfield-Latimer/Gow-Devaux models.
                 
Edit History: Created by Nilapratim Sengupta in July-August 2023.
"""
def proc_setBiophysics_pyramidalAxon(pyramidalCell, list_extracellularParameters_axon):
    
    for sec in pyramidalCell.pyramidalAxon.TotalAxon:
        sec.Ra = 100  # (ohm.cm) Gets updated later
        sec.cm = 0.9  # (uF/cm^2) Gets updated later
        sec.insert('pas')
        sec.g_pas = 0.001  # Gets updated later
        sec.e_pas = -70

    for sec in pyramidalCell.pyramidalAxon.TotalAxon:
        sec.insert('naf')
        sec.insert('kdr')
        sec.insert('ka')
        sec.insert('k2')

        sec.gbar_naf = 400e-3  # Gets updated later
        sec.gbar_kdr = 400e-3  # Gets updated later
        sec.gbar_ka = 2e-3  # Gets updated later
        sec.gbar_k2 = 0.1e-3  # Gets updated later

    for sec in pyramidalCell.pyramidalAxon.Nodes_axon:
        sec.insert('extracellular')
        sec.xraxial[0] = list_extracellularParameters_axon[0][0]  # (Mohm/cm)
        sec.xraxial[1] = list_extracellularParameters_axon[1][0]  # (Mohm/cm)
        sec.xg[0] = list_extracellularParameters_axon[2][0]  # (S/cm^2)
        sec.xg[1] = list_extracellularParameters_axon[3][0]  # (S/cm^2)
        sec.xc[0] = list_extracellularParameters_axon[4][0]  # (uF/cm^2)
        sec.xc[1] = list_extracellularParameters_axon[5][0]  # (uF/cm^2)
        sec.e_extracellular = list_extracellularParameters_axon[6][0]  # (mV)

    for sec in pyramidalCell.pyramidalAxon.ParanodeOnes_axon:
        sec.insert('extracellular')
        sec.xraxial[0] = list_extracellularParameters_axon[0][1]  # (Mohm/cm)
        sec.xraxial[1] = list_extracellularParameters_axon[1][1]  # (Mohm/cm)
        sec.xg[0] = list_extracellularParameters_axon[2][1]  # (S/cm^2)
        sec.xg[1] = list_extracellularParameters_axon[3][1]  # (S/cm^2)
        sec.xc[0] = list_extracellularParameters_axon[4][1]  # (uF/cm^2)
        sec.xc[1] = list_extracellularParameters_axon[5][1]  # (uF/cm^2)
        sec.e_extracellular = list_extracellularParameters_axon[6][1]  # (mV)

    for sec in pyramidalCell.pyramidalAxon.ParanodeTwos_axon:
        sec.insert('extracellular')
        sec.xraxial[0] = list_extracellularParameters_axon[0][2]  # (Mohm/cm)
        sec.xraxial[1] = list_extracellularParameters_axon[1][2]  # (Mohm/cm)
        sec.xg[0] = list_extracellularParameters_axon[2][2]  # (S/cm^2)
        sec.xg[1] = list_extracellularParameters_axon[3][2]  # (S/cm^2)
        sec.xc[0] = list_extracellularParameters_axon[4][2]  # (uF/cm^2)
        sec.xc[1] = list_extracellularParameters_axon[5][2]  # (uF/cm^2)
        sec.e_extracellular = list_extracellularParameters_axon[6][2]  # (mV)
        
    for sec in pyramidalCell.pyramidalAxon.ParanodeThrees_axon:
        sec.insert('extracellular')
        sec.xraxial[0] = list_extracellularParameters_axon[0][3]  # (Mohm/cm)
        sec.xraxial[1] = list_extracellularParameters_axon[1][3]  # (Mohm/cm)
        sec.xg[0] = list_extracellularParameters_axon[2][3]  # (S/cm^2)
        sec.xg[1] = list_extracellularParameters_axon[3][3]  # (S/cm^2)
        sec.xc[0] = list_extracellularParameters_axon[4][3]  # (uF/cm^2)
        sec.xc[1] = list_extracellularParameters_axon[5][3]  # (uF/cm^2)
        sec.e_extracellular = list_extracellularParameters_axon[6][3]  # (mV)     
        
    for sec in pyramidalCell.pyramidalAxon.ParanodeFours_axon:
        sec.insert('extracellular')
        sec.xraxial[0] = list_extracellularParameters_axon[0][4]  # (Mohm/cm)
        sec.xraxial[1] = list_extracellularParameters_axon[1][4]  # (Mohm/cm)
        sec.xg[0] = list_extracellularParameters_axon[2][4]  # (S/cm^2)
        sec.xg[1] = list_extracellularParameters_axon[3][4]  # (S/cm^2)
        sec.xc[0] = list_extracellularParameters_axon[4][4]  # (uF/cm^2)
        sec.xc[1] = list_extracellularParameters_axon[5][4]  # (uF/cm^2)
        sec.e_extracellular = list_extracellularParameters_axon[6][4]  # (mV)   
        
    for sec in pyramidalCell.pyramidalAxon.Juxtaparanodes_axon:
        sec.insert('extracellular')
        sec.xraxial[0] = list_extracellularParameters_axon[0][5]  # (Mohm/cm)
        sec.xraxial[1] = list_extracellularParameters_axon[1][5]  # (Mohm/cm)
        sec.xg[0] = list_extracellularParameters_axon[2][5]  # (S/cm^2)
        sec.xg[1] = list_extracellularParameters_axon[3][5]  # (S/cm^2)
        sec.xc[0] = list_extracellularParameters_axon[4][5]  # (uF/cm^2)
        sec.xc[1] = list_extracellularParameters_axon[5][5]  # (uF/cm^2)
        sec.e_extracellular = list_extracellularParameters_axon[6][5]  # (mV)
        
    for sec in pyramidalCell.pyramidalAxon.Internodes_axon:
        sec.insert('extracellular')
        sec.xraxial[0] = list_extracellularParameters_axon[0][6]  # (Mohm/cm)
        sec.xraxial[1] = list_extracellularParameters_axon[1][6]  # (Mohm/cm)
        sec.xg[0] = list_extracellularParameters_axon[2][6]  # (S/cm^2)
        sec.xg[1] = list_extracellularParameters_axon[3][6]  # (S/cm^2)
        sec.xc[0] = list_extracellularParameters_axon[4][6]  # (uF/cm^2)
        sec.xc[1] = list_extracellularParameters_axon[5][6]  # (uF/cm^2)
        sec.e_extracellular = list_extracellularParameters_axon[6][6]  # (mV)
        
# End of proc_setBiophysics_pyramidalAxon()





"""
Name: proc_setBiophysics_pyramidalCollateral()
Type: Procedure

Input Parameter(s): 1. Object - pyramidalCell
                    2. Index of collateral
                    3. List of lists (of parameter values for extracellular implementation) - [[xraxial0],[xraxial1],[xg0],[xg1],[xc0],[xc1],[e]]
Return Parameter: None

Description: Procedure to set active and passive membrane mechanism properties for different sections of the collateral(s) of the pyramidal cell.
             The corresponding values for the collateral sections are calculated based on the Scurfield-Latimer/Gow-Devaux models.
                 
Edit History: Created by Nilapratim Sengupta in July-August 2023.
              Edited by Nilapratim Sengupta in September 2024 to ensure collaterals can be modified individually.
"""
def proc_setBiophysics_pyramidalCollateral(pyramidalCell, index_collateral, list_extracellularParameters_collateral):
    
    for sec in pyramidalCell.pyramidalCollateral[index_collateral].TotalCollateral:
        sec.Ra = 100  # (ohm.cm) Gets updated later
        sec.cm = 0.9  # (uF/cm^2) Gets updated later
        sec.insert('pas')
        sec.g_pas = 0.001  # Gets updated later
        sec.e_pas = -70

    for sec in pyramidalCell.pyramidalCollateral[index_collateral].TotalCollateral:
        sec.insert('naf')
        sec.insert('kdr')
        sec.insert('ka')
        sec.insert('k2')

        sec.gbar_naf = 400e-3  # Gets updated later
        sec.gbar_kdr = 400e-3  # Gets updated later
        sec.gbar_ka = 2e-3  # Gets updated later
        sec.gbar_k2 = 0.1e-3  # Gets updated later

    for sec in pyramidalCell.pyramidalCollateral[index_collateral].Nodes_collateral:
        sec.insert('extracellular')
        sec.xraxial[0] = list_extracellularParameters_collateral[0][0]  # (Mohm/cm)
        sec.xraxial[1] = list_extracellularParameters_collateral[1][0]  # (Mohm/cm)
        sec.xg[0] = list_extracellularParameters_collateral[2][0]  # (S/cm^2)
        sec.xg[1] = list_extracellularParameters_collateral[3][0]  # (S/cm^2)
        sec.xc[0] = list_extracellularParameters_collateral[4][0]  # (uF/cm^2)
        sec.xc[1] = list_extracellularParameters_collateral[5][0]  # (uF/cm^2)
        sec.e_extracellular = list_extracellularParameters_collateral[6][0]  # (mV)

    for sec in pyramidalCell.pyramidalCollateral[index_collateral].ParanodeOnes_collateral:
        sec.insert('extracellular')
        sec.xraxial[0] = list_extracellularParameters_collateral[0][1]  # (Mohm/cm)
        sec.xraxial[1] = list_extracellularParameters_collateral[1][1]  # (Mohm/cm)
        sec.xg[0] = list_extracellularParameters_collateral[2][1]  # (S/cm^2)
        sec.xg[1] = list_extracellularParameters_collateral[3][1]  # (S/cm^2)
        sec.xc[0] = list_extracellularParameters_collateral[4][1]  # (uF/cm^2)
        sec.xc[1] = list_extracellularParameters_collateral[5][1]  # (uF/cm^2)
        sec.e_extracellular = list_extracellularParameters_collateral[6][1]  # (mV)

    for sec in pyramidalCell.pyramidalCollateral[index_collateral].ParanodeTwos_collateral:
        sec.insert('extracellular')
        sec.xraxial[0] = list_extracellularParameters_collateral[0][2]  # (Mohm/cm)
        sec.xraxial[1] = list_extracellularParameters_collateral[1][2]  # (Mohm/cm)
        sec.xg[0] = list_extracellularParameters_collateral[2][2]  # (S/cm^2)
        sec.xg[1] = list_extracellularParameters_collateral[3][2]  # (S/cm^2)
        sec.xc[0] = list_extracellularParameters_collateral[4][2]  # (uF/cm^2)
        sec.xc[1] = list_extracellularParameters_collateral[5][2]  # (uF/cm^2)
        sec.e_extracellular = list_extracellularParameters_collateral[6][2]  # (mV)
        
    for sec in pyramidalCell.pyramidalCollateral[index_collateral].ParanodeThrees_collateral:
        sec.insert('extracellular')
        sec.xraxial[0] = list_extracellularParameters_collateral[0][3]  # (Mohm/cm)
        sec.xraxial[1] = list_extracellularParameters_collateral[1][3]  # (Mohm/cm)
        sec.xg[0] = list_extracellularParameters_collateral[2][3]  # (S/cm^2)
        sec.xg[1] = list_extracellularParameters_collateral[3][3]  # (S/cm^2)
        sec.xc[0] = list_extracellularParameters_collateral[4][3]  # (uF/cm^2)
        sec.xc[1] = list_extracellularParameters_collateral[5][3]  # (uF/cm^2)
        sec.e_extracellular = list_extracellularParameters_collateral[6][3]  # (mV)     
        
    for sec in pyramidalCell.pyramidalCollateral[index_collateral].ParanodeFours_collateral:
        sec.insert('extracellular')
        sec.xraxial[0] = list_extracellularParameters_collateral[0][4]  # (Mohm/cm)
        sec.xraxial[1] = list_extracellularParameters_collateral[1][4]  # (Mohm/cm)
        sec.xg[0] = list_extracellularParameters_collateral[2][4]  # (S/cm^2)
        sec.xg[1] = list_extracellularParameters_collateral[3][4]  # (S/cm^2)
        sec.xc[0] = list_extracellularParameters_collateral[4][4]  # (uF/cm^2)
        sec.xc[1] = list_extracellularParameters_collateral[5][4]  # (uF/cm^2)
        sec.e_extracellular = list_extracellularParameters_collateral[6][4]  # (mV)   
        
    for sec in pyramidalCell.pyramidalCollateral[index_collateral].Juxtaparanodes_collateral:
        sec.insert('extracellular')
        sec.xraxial[0] = list_extracellularParameters_collateral[0][5]  # (Mohm/cm)
        sec.xraxial[1] = list_extracellularParameters_collateral[1][5]  # (Mohm/cm)
        sec.xg[0] = list_extracellularParameters_collateral[2][5]  # (S/cm^2)
        sec.xg[1] = list_extracellularParameters_collateral[3][5]  # (S/cm^2)
        sec.xc[0] = list_extracellularParameters_collateral[4][5]  # (uF/cm^2)
        sec.xc[1] = list_extracellularParameters_collateral[5][5]  # (uF/cm^2)
        sec.e_extracellular = list_extracellularParameters_collateral[6][5]  # (mV)
        
    for sec in pyramidalCell.pyramidalCollateral[index_collateral].Internodes_collateral:
        sec.insert('extracellular')
        sec.xraxial[0] = list_extracellularParameters_collateral[0][6]  # (Mohm/cm)
        sec.xraxial[1] = list_extracellularParameters_collateral[1][6]  # (Mohm/cm)
        sec.xg[0] = list_extracellularParameters_collateral[2][6]  # (S/cm^2)
        sec.xg[1] = list_extracellularParameters_collateral[3][6]  # (S/cm^2)
        sec.xc[0] = list_extracellularParameters_collateral[4][6]  # (uF/cm^2)
        sec.xc[1] = list_extracellularParameters_collateral[5][6]  # (uF/cm^2)
        sec.e_extracellular = list_extracellularParameters_collateral[6][6]  # (mV)
        
# End of proc_setBiophysics_pyramidalCollateral()





"""
Name: proc_applySpineCorrection()
Type: Procedure

Input Parameter(s): 1. Object - pyramidalCell
                    2. Method of spine correction (either 1 or 2)
Return Parameter: None

Description: Procedure to make necessary adjustments to the dendritic properties to account for spines.
             The corresponding values for the dendrites are from the Rumbell/Traub models.
                 
Edit History: Created by Nilapratim Sengupta in July-August 2023.
"""
def proc_applySpineCorrection(pyramidalCell, method_spineCorrection):
    
    # Adjusting dendritic properties to account for spines
    for sec in pyramidalCell.pyramidalSomaDendrites.Dendrites:
        
        if (method_spineCorrection == 1):
            sec.L *= 2
            sec.Ra /= 2 # Check later
        elif (method_spineCorrection == 2):
            sec.g_pas *= 2
            sec.cm *= 2
            sec.gbar_naf *= 2
            sec.gbar_nap *= 2
            sec.gbar_kdr *= 2
            sec.gbar_ka *= 2   
            sec.gbar_k2 *= 2
            sec.gbar_km *= 2
            sec.gbar_ar *= 2
            sec.gbar_kc *= 2
            sec.gbar_kahp *= 2
            sec.gbar_cat *= 2
            sec.gbar_cal *= 2
            sec.phi_cad /= 2
            
    
# End of proc_applySpineCorrection()





"""
Name: proc_applyPerturbation_pyramidalAxon()
Type: Procedure

Input Parameter(s): 1. Object - pyramidal axon
                    2. Index of affected segment (internode)
                    3. Diameter of the axon (um)
                    4. Length of each node um)
                    5. Length of each myelinated segment (um)
                    6. Number of myelin lamellae over each segment - perturbed
                    7. Thickness of each myelin lamella (um)
                    8. Perturbation code (0: unperturbed, 1: demyelinated, >=2: remyelinated)
                    9. Paranode-scaling code (0: fixed, 1: scaled)
Return Parameter: None

Description: Procedure to update sections of perturbed segment in the myelinated axon.
             This procedure is based on calculations listed in an Excel sheet, prepared by Gow-Devaux.


Edit History: Created by Nilapratim Sengupta in August 2023.
"""
def proc_applyPerturbation_pyramidalAxon(pyramidalObject, index_affectedSegment, diameter_axon, length_node, length_segment, count_lamellae, thickness_lamella, code_perturbation, code_paranodeScaling):

    # Lengths - compartmental and total
    [compartmentalLength_paranode, compartmentalLength_juxtaparanode, totalLength_internode] = func_calculateSectionLengths(diameter_axon, length_node, length_segment, code_perturbation, code_paranodeScaling)
     
    # Calculating values of extracellular parameters for the perturbed segment
    list_extracellularParameters_axon = func_calculateExtracellularParameters(diameter_axon, length_node, length_segment, count_lamellae, thickness_lamella, code_perturbation, code_paranodeScaling)

    # Applying updated values for parameters to implement the perturbation
    for loopCounter in range(2*index_affectedSegment, 2*index_affectedSegment+2): # Note: In Python loopCounter will go until (2*index_affectedSegment+2)-1

        for index, sec in enumerate(pyramidalObject.ParanodeOnes_axon):
            if index == loopCounter:
                sec.L = compartmentalLength_paranode
                sec.xraxial[0] = list_extracellularParameters_axon[0][1]  # (Mohm/cm)
                sec.xraxial[1] = list_extracellularParameters_axon[1][1]  # (Mohm/cm)
                sec.xg[0] = list_extracellularParameters_axon[2][1]  # (S/cm^2)
                sec.xg[1] = list_extracellularParameters_axon[3][1]  # (S/cm^2)
                sec.xc[0] = list_extracellularParameters_axon[4][1]  # (uF/cm^2)
                sec.xc[1] = list_extracellularParameters_axon[5][1]  # (uF/cm^2)
                sec.e_extracellular = list_extracellularParameters_axon[6][1]  # (mV)

        for index, sec in enumerate(pyramidalObject.ParanodeTwos_axon):
            if index == loopCounter:
                sec.L = compartmentalLength_paranode
                sec.xraxial[0] = list_extracellularParameters_axon[0][2]  # (Mohm/cm)
                sec.xraxial[1] = list_extracellularParameters_axon[1][2]  # (Mohm/cm)
                sec.xg[0] = list_extracellularParameters_axon[2][2]  # (S/cm^2)
                sec.xg[1] = list_extracellularParameters_axon[3][2]  # (S/cm^2)
                sec.xc[0] = list_extracellularParameters_axon[4][2]  # (uF/cm^2)
                sec.xc[1] = list_extracellularParameters_axon[5][2]  # (uF/cm^2)
                sec.e_extracellular = list_extracellularParameters_axon[6][2]  # (mV)
            
        for index, sec in enumerate(pyramidalObject.ParanodeThrees_axon):
            if index == loopCounter:
                sec.L = compartmentalLength_paranode
                sec.xraxial[0] = list_extracellularParameters_axon[0][3]  # (Mohm/cm)
                sec.xraxial[1] = list_extracellularParameters_axon[1][3]  # (Mohm/cm)
                sec.xg[0] = list_extracellularParameters_axon[2][3]  # (S/cm^2)
                sec.xg[1] = list_extracellularParameters_axon[3][3]  # (S/cm^2)
                sec.xc[0] = list_extracellularParameters_axon[4][3]  # (uF/cm^2)
                sec.xc[1] = list_extracellularParameters_axon[5][3]  # (uF/cm^2)
                sec.e_extracellular = list_extracellularParameters_axon[6][3]  # (mV)     
            
        for index, sec in enumerate(pyramidalObject.ParanodeFours_axon):
            if index == loopCounter:
                sec.L = compartmentalLength_paranode
                sec.xraxial[0] = list_extracellularParameters_axon[0][4]  # (Mohm/cm)
                sec.xraxial[1] = list_extracellularParameters_axon[1][4]  # (Mohm/cm)
                sec.xg[0] = list_extracellularParameters_axon[2][4]  # (S/cm^2)
                sec.xg[1] = list_extracellularParameters_axon[3][4]  # (S/cm^2)
                sec.xc[0] = list_extracellularParameters_axon[4][4]  # (uF/cm^2)
                sec.xc[1] = list_extracellularParameters_axon[5][4]  # (uF/cm^2)
                sec.e_extracellular = list_extracellularParameters_axon[6][4]  # (mV)   
            
        for index, sec in enumerate(pyramidalObject.Juxtaparanodes_axon):
            if index == loopCounter:
                sec.L = compartmentalLength_juxtaparanode
                sec.xraxial[0] = list_extracellularParameters_axon[0][5]  # (Mohm/cm)
                sec.xraxial[1] = list_extracellularParameters_axon[1][5]  # (Mohm/cm)
                sec.xg[0] = list_extracellularParameters_axon[2][5]  # (S/cm^2)
                sec.xg[1] = list_extracellularParameters_axon[3][5]  # (S/cm^2)
                sec.xc[0] = list_extracellularParameters_axon[4][5]  # (uF/cm^2)
                sec.xc[1] = list_extracellularParameters_axon[5][5]  # (uF/cm^2)
                sec.e_extracellular = list_extracellularParameters_axon[6][5]  # (mV)
        
    for index, sec in enumerate(pyramidalObject.Internodes_axon):
        if index == index_affectedSegment:
            sec.L = totalLength_internode
            sec.xraxial[0] = list_extracellularParameters_axon[0][6]  # (Mohm/cm)
            sec.xraxial[1] = list_extracellularParameters_axon[1][6]  # (Mohm/cm)
            sec.xg[0] = list_extracellularParameters_axon[2][6]  # (S/cm^2)
            sec.xg[1] = list_extracellularParameters_axon[3][6]  # (S/cm^2)
            sec.xc[0] = list_extracellularParameters_axon[4][6]  # (uF/cm^2)
            sec.xc[1] = list_extracellularParameters_axon[5][6]  # (uF/cm^2)
            sec.e_extracellular = list_extracellularParameters_axon[6][6]  # (mV)
    
# End of proc_applyPerturbation_pyramidalAxon()





"""
Name: proc_applyPerturbation_pyramidalCollateral()
Type: Procedure

Input Parameter(s): 1. Object - pyramidal collateral
                    2. Index of affected segment (internode)
                    3. Diameter of the collateral (um)
                    4. Length of each node um)
                    5. Length of each myelinated segment (um)
                    6. Number of myelin lamellae over each segment - perturbed
                    7. Thickness of each myelin lamella (um)
                    8. Perturbation code (0: unperturbed, 1: demyelinated, >=2: remyelinated)
                    9. Paranode-scaling code (0: fixed, 1: scaled)
Return Parameter: None

Description: Procedure to update sections of perturbed segment in a myelinated collateral.
             This procedure is based on calculations listed in an Excel sheet, prepared by Gow-Devaux.


Edit History: Created by Nilapratim Sengupta in August 2023.
              Edited by Nilapratim Sengupta in January 2024 to update collateral diameter.
"""
def proc_applyPerturbation_pyramidalCollateral(pyramidalObject, index_affectedSegment, diameter_collateral, length_node, length_segment, count_lamellae, thickness_lamella, code_perturbation, code_paranodeScaling):

    # Lengths - compartmental and total
    [compartmentalLength_paranode, compartmentalLength_juxtaparanode, totalLength_internode] = func_calculateSectionLengths(diameter_collateral, length_node, length_segment, code_perturbation, code_paranodeScaling)
     
    # Calculating values of extracellular parameters for the perturbed segment
    list_extracellularParameters_collateral = func_calculateExtracellularParameters(diameter_collateral, length_node, length_segment, count_lamellae, thickness_lamella, code_perturbation, code_paranodeScaling)
    
    # Applying updated values for parameters to implement the perturbation
    for loopCounter in range(2*index_affectedSegment, 2*index_affectedSegment+2): # Note: In Python loopCounter will go until (2*index_affectedSegment+2)-1

        for index, sec in enumerate(pyramidalObject.ParanodeOnes_collateral):
            if index == loopCounter:
                sec.L = compartmentalLength_paranode
                sec.xraxial[0] = list_extracellularParameters_collateral[0][1]  # (Mohm/cm)
                sec.xraxial[1] = list_extracellularParameters_collateral[1][1]  # (Mohm/cm)
                sec.xg[0] = list_extracellularParameters_collateral[2][1]  # (S/cm^2)
                sec.xg[1] = list_extracellularParameters_collateral[3][1]  # (S/cm^2)
                sec.xc[0] = list_extracellularParameters_collateral[4][1]  # (uF/cm^2)
                sec.xc[1] = list_extracellularParameters_collateral[5][1]  # (uF/cm^2)
                sec.e_extracellular = list_extracellularParameters_collateral[6][1]  # (mV)

        for index, sec in enumerate(pyramidalObject.ParanodeTwos_collateral):
            if index == loopCounter:
                sec.L = compartmentalLength_paranode
                sec.xraxial[0] = list_extracellularParameters_collateral[0][2]  # (Mohm/cm)
                sec.xraxial[1] = list_extracellularParameters_collateral[1][2]  # (Mohm/cm)
                sec.xg[0] = list_extracellularParameters_collateral[2][2]  # (S/cm^2)
                sec.xg[1] = list_extracellularParameters_collateral[3][2]  # (S/cm^2)
                sec.xc[0] = list_extracellularParameters_collateral[4][2]  # (uF/cm^2)
                sec.xc[1] = list_extracellularParameters_collateral[5][2]  # (uF/cm^2)
                sec.e_extracellular = list_extracellularParameters_collateral[6][2]  # (mV)
            
        for index, sec in enumerate(pyramidalObject.ParanodeThrees_collateral):
            if index == loopCounter:
                sec.L = compartmentalLength_paranode
                sec.xraxial[0] = list_extracellularParameters_collateral[0][3]  # (Mohm/cm)
                sec.xraxial[1] = list_extracellularParameters_collateral[1][3]  # (Mohm/cm)
                sec.xg[0] = list_extracellularParameters_collateral[2][3]  # (S/cm^2)
                sec.xg[1] = list_extracellularParameters_collateral[3][3]  # (S/cm^2)
                sec.xc[0] = list_extracellularParameters_collateral[4][3]  # (uF/cm^2)
                sec.xc[1] = list_extracellularParameters_collateral[5][3]  # (uF/cm^2)
                sec.e_extracellular = list_extracellularParameters_collateral[6][3]  # (mV)     
            
        for index, sec in enumerate(pyramidalObject.ParanodeFours_collateral):
            if index == loopCounter:
                sec.L = compartmentalLength_paranode
                sec.xraxial[0] = list_extracellularParameters_collateral[0][4]  # (Mohm/cm)
                sec.xraxial[1] = list_extracellularParameters_collateral[1][4]  # (Mohm/cm)
                sec.xg[0] = list_extracellularParameters_collateral[2][4]  # (S/cm^2)
                sec.xg[1] = list_extracellularParameters_collateral[3][4]  # (S/cm^2)
                sec.xc[0] = list_extracellularParameters_collateral[4][4]  # (uF/cm^2)
                sec.xc[1] = list_extracellularParameters_collateral[5][4]  # (uF/cm^2)
                sec.e_extracellular = list_extracellularParameters_collateral[6][4]  # (mV)   
            
        for index, sec in enumerate(pyramidalObject.Juxtaparanodes_collateral):
            if index == loopCounter:
                sec.L = compartmentalLength_juxtaparanode
                sec.xraxial[0] = list_extracellularParameters_collateral[0][5]  # (Mohm/cm)
                sec.xraxial[1] = list_extracellularParameters_collateral[1][5]  # (Mohm/cm)
                sec.xg[0] = list_extracellularParameters_collateral[2][5]  # (S/cm^2)
                sec.xg[1] = list_extracellularParameters_collateral[3][5]  # (S/cm^2)
                sec.xc[0] = list_extracellularParameters_collateral[4][5]  # (uF/cm^2)
                sec.xc[1] = list_extracellularParameters_collateral[5][5]  # (uF/cm^2)
                sec.e_extracellular = list_extracellularParameters_collateral[6][5]  # (mV)
        
    for index, sec in enumerate(pyramidalObject.Internodes_collateral):
        if index == index_affectedSegment:
            sec.L = totalLength_internode
            sec.xraxial[0] = list_extracellularParameters_collateral[0][6]  # (Mohm/cm)
            sec.xraxial[1] = list_extracellularParameters_collateral[1][6]  # (Mohm/cm)
            sec.xg[0] = list_extracellularParameters_collateral[2][6]  # (S/cm^2)
            sec.xg[1] = list_extracellularParameters_collateral[3][6]  # (S/cm^2)
            sec.xc[0] = list_extracellularParameters_collateral[4][6]  # (uF/cm^2)
            sec.xc[1] = list_extracellularParameters_collateral[5][6]  # (uF/cm^2)
            sec.e_extracellular = list_extracellularParameters_collateral[6][6]  # (mV)
    
# End of proc_applyPerturbation_pyramidalCollateral()





"""
Name: func_setupIClamp()
Type: Function

Input Parameter(s): 1. Object - pyramidalCell
                    2. Start time of current clamp
                    3. Stop time of current clamp
Return Parameter: 1. List of objects - holdingStimulus, currentStimulus

Description: Function to create holding stimulus and current stimulus, applied to the center of the soma.
                 
Edit History: Created by Nilapratim Sengupta in July-August 2023.
"""
def func_setupIClamp(pyramidalCell, time_currentStart, time_currentStop):
    
    # Holding phase
    holdingStimulus = h.IClamp(0.5, sec=pyramidalCell.pyramidalSomaDendrites.comp[1])
    holdingStimulus.dur = 2015
    holdingStimulus.delay = -2000
    holdingStimulus.amp = 0

    # Stimulation phase
    currentStimulus = h.IClamp(0.5, sec=pyramidalCell.pyramidalSomaDendrites.comp[1])
    currentStimulus.dur = time_currentStop - time_currentStart  
    currentStimulus.delay = time_currentStart
    currentStimulus.amp = -0.15 # Gets modified later

    return [holdingStimulus, currentStimulus]
    
# End of func_setupIClamp()





"""
Name: func_setupVClamp()
Type: Function

Input Parameter(s): 1. Object - pyramidalCell
                    2. Initial membrane potential
Return Parameter: 1. Object - voltageStimulus

Description: Function to create instance of a voltage clamp (defined in vsource.mod), applied to the center of the soma.
                 
Edit History: Created by Nilapratim Sengupta in August 2023.
"""
def func_setupVClamp(pyramidalCell, membranePotential_initial):
    
    # Voltage clamp
    voltageStimulus = h.Vsource(0.5, sec=pyramidalCell.pyramidalSomaDendrites.comp[1])
    voltageStimulus.rs = 0.01
    voltageStimulus.toff = 0   
    voltageStimulus.amp = membranePotential_initial

    return voltageStimulus
    
# End of func_setupVClamp()





"""
Name: func_loadChannelParameters()
Type: Function

Input Parameter(s): 1. Type of cell (1: default, 2: young; 3: aged)
                    2. Parameter index (Either 0, 1 or 2)
Return Parameter: 1. List of ion channel parameters 

Description: Function to load ion channel parameters in accordance with Rumbell/Traub models.
                 
Edit History: Created by Nilapratim Sengupta in July-August 2023.
"""
def func_loadChannelParameters(type_cell, index_parameter):
    
    # Creating and loading lists with parameter values, based on type of cell
    if (type_cell >= 2):
        list_parametersZero = [0.042039, 0.31611, 0.0047585, 0.0106580, 0.0037284, 0.00015574, 9.0673e-05, 830.25, 0.0097944, 1.7771, 8.3288, 18.631, 4.1197, 9.4392, 6.0591, 13.646, 1.1585, -15.304, -14.901, -29.992, 0.43005, -36.416, -29.279]
        list_parametersOne = [0.042802, 0.38121, 0.0074711, 0.0065173, 0.0036778, 0.00011990, 0.00011110, 736.04, 0.0099149, 2.0840, 8.2943, 18.156, 3.0562, 12.087, 15.652, 12.510, 1.4208, -15.248, -16.265, -34.251, 3.00000, -34.010, -31.586]
        list_parametersTwo = [0.041768, 0.14555, 0.0082385, 0.0102910, 0.0050655, 6.6648e-05, 8.8751e-05, 298.08, 0.0097699, 1.6815, 7.2567, 18.566, 2.5945, 16.806, 11.124, 14.507, 1.7094, -15.111, -11.650, -35.180, -21.203, -27.749, -21.756]
    elif (type_cell == 3):
        list_parametersZero = [0.088676, 0.282090, 0.0054771, 0.019922, 0.00065710, 6.0722e-06, 0.0009996499999999999, 47479, 0.0020877, 6.2912, 2.3057, 3.2847, 2.3130, 1.3855, 4.7349, 11.179, 1.9527, -5.4739, -36.983, -4.8420, -15.314, 10.5170, -10.809]
        list_parametersOne = [0.021656, 0.022436, 0.0127170, 0.016599, 0.00093677, 0.00029288, 0.0005505300000000000, 34627, 0.0097155, 4.0338, 3.2426, 10.043, 11.298, 14.579, 11.132, 16.244, 1.0253, -8.8705, -10.170, -21.948, -19.727, -35.308, -32.771]
        list_parametersTwo = [0.088802, 0.217980, 0.0049884, 0.019907, 0.00067711, 4.5122e-06, 0.0009936200000000000, 51753, 0.0026811, 6.2350, 2.3197, 9.0452, 2.2796, 1.4942, 6.0396, 8.3907, 1.9444, -4.9373, -34.492, -4.3734, -14.967, 10.1060, -10.889]

    # Returning desired list of ion channel parameters based on parameter-index
    match index_parameter:
        case 0:
            return list_parametersZero
        case 1:
            return list_parametersOne
        case 2:
            return list_parametersTwo
   
# End of func_loadChannelParameters()





"""
Name: proc_setChannelConductances_pyramidalSomaDendrites()
Type: Procedure

Input Parameter(s): 1. Object - pyramidalCell
                    2. Type of cell (1: default, 2: young; 3: aged)
                    3. Parameter index (Either 0, 1 or 2)
Return Parameter: None

Description: Procedure to set ion channel conductances in accordance with Rumbell/Traub models for the soma and dendrites.
                 
Edit History: Created by Nilapratim Sengupta in July-August 2023.
"""
def proc_setChannelConductances_pyramidalSomaDendrites(pyramidalCell, type_cell, index_parameter):
    
    # Dummy conductance parameters for ion channels
    dphi_cad   = 52 / 2e-3
    dbeta_cad  = 1 / 100	# In the paper beta = 50 [ms]

    dgnaf	= 0.2501575996 # Tuned NernstTest12 // 0.312 // tuned 84 // 150e-3 * 1.25 // default
    dgkdr  	= 0.0513842372 # Tuned NernstTest12 // 1.554 // tuned 84 // 125e-3
    dgka   	= 0
    dgk2   	= 0
    dgkm	= 0.0016943554 # Tuned NernstTest12 // 0.00138 // tuned 84 // 2.5 * 1.5e-3 * 2
    dgkc	= 0
    dgkahp	= 0

    dgcad 	= 0
    dgnap 	= 0
    dgcat 	= 0
    dgcal 	= 0
    
    # Dummy parameters for passive channels
    if (type_cell == 2):
        dgar = 2.83162e-05 # 1e-10
        dg_pas = 3.95094e-05 # 9.0916e-05
        de_pas = -67.92914 # -69.9598
        dcm = 1.09243 # 1.0 // (manual) // 0.9 // (default)
    elif (type_cell == 3):
        dgar = 3.61324e-05 # 5.0560e-05
        dg_pas = 1.71951e-05 # 1.7479e-05
        de_pas = -76.7721 # -79.9509
        dcm = 0.70449 # 1.0 // (manual) // 0.9 // (default)

    dRa	  	= 150	# (manual) // 250 // (default)
    dena	= 65.204611 # Using nernst( 115 , 5 , 1 )  // 50 // (default)
    dek		= -75.497 # Using nernst( 115 , 5 , 1 ) // -95 // (default)
    
    # Acquiring the list of ion channel parameter values based on type of cell and parameter index
    list_channelParameters = func_loadChannelParameters(type_cell, index_parameter)
    
    # Assigning values to designated parameters
    dgnaf = list_channelParameters[0]
    dgkdr = list_channelParameters[1]
    dgkm = list_channelParameters[2]
    dgka = list_channelParameters[3]
    dgkahp = list_channelParameters[4]
    dgnap = list_channelParameters[5]
    dgcal = list_channelParameters[6]
    dphi_cad = list_channelParameters[7]
    dbeta_cad = list_channelParameters[8]
    dtmmnaf = list_channelParameters[9]
    dtmhnaf = list_channelParameters[10]
    dtmkdr = list_channelParameters[11]
    dtmkm = list_channelParameters[12]
    dtmmka = list_channelParameters[13]
    dtmhka = list_channelParameters[14]
    dtmnap = list_channelParameters[15]
    dtmcal = list_channelParameters[16]
    dvsnaf = list_channelParameters[17]
    dvskdr = list_channelParameters[18]
    dvskm = list_channelParameters[19]
    dvska = list_channelParameters[20]
    dvsnap = list_channelParameters[21]
    dvscal = list_channelParameters[22]

 
    # Assigning values to ion channel conductances for the soma and dendrites
    for sec in pyramidalCell.pyramidalSomaDendrites.Soma:
        sec.phi_cad = dphi_cad
        sec.beta_cad = dbeta_cad

        sec.gbar_naf = dgnaf
        sec.gbar_kdr = dgkdr
        sec.gbar_ka = dgka
        sec.gbar_k2 = dgk2
        sec.gbar_ar = dgar
        sec.gbar_kc = 1.6 * dgkc  # dkc * dgkc  # In the paper 'dkc * 12e-3'
        sec.gbar_kahp = dgkahp
        sec.gbar_km = dgkm

        # sec.gbar_cad = dgcad # Check later
        sec.gbar_nap = dgnap
        sec.gbar_cat = dgcat
        sec.gbar_cal = dgcal

        sec.g_pas = dg_pas
        sec.e_pas = de_pas
        sec.Ra = dRa
        sec.cm = dcm
        sec.ena = dena
        sec.ek = dek

    for sec in pyramidalCell.pyramidalSomaDendrites.Dendrites:
        sec.phi_cad = dphi_cad
        sec.beta_cad = dbeta_cad * 5

        sec.gbar_naf = dgnaf * 0.03333
        sec.gbar_kdr = 0
        sec.gbar_ka = dgka * 0.06667
        sec.gbar_k2 = dgk2
        sec.gbar_ar = dgar
        sec.gbar_kc = 0
        sec.gbar_kahp = dgkahp
        sec.gbar_km = dgkm
     
        # sec.gbar_cad = dgcad # Check later
        sec.gbar_nap = dgnap
        sec.gbar_cat = dgcat
        sec.gbar_cal = dgcal

        sec.g_pas = dg_pas
        sec.e_pas = de_pas
        sec.Ra = dRa
        sec.cm = dcm
        sec.ena = dena
        sec.ek = dek

    for sec in pyramidalCell.pyramidalSomaDendrites.Proximal:
        sec.gbar_naf = dgnaf * 0.5
        sec.gbar_kdr = dgkdr * 0.75

    pyramidalCell.pyramidalSomaDendrites.comp[38].gbar_ka = dgka
    pyramidalCell.pyramidalSomaDendrites.comp[38].gbar_naf = dgnaf * 0.66667
    pyramidalCell.pyramidalSomaDendrites.comp[38].gbar_kdr = dgkdr * 0.75        
    

# End of proc_setChannelConductances_pyramidalSomaDendrites()





"""
Name: proc_setChannelConductances_pyramidalAxon()
Type: Procedure

Input Parameter(s): 1. Object - pyramidalCell
                    2. Type of cell (1: default, 2: young; 3: aged)
                    3. Parameter index (Either 0, 1 or 2)
Return Parameter: None

Description: Procedure to set ion channel conductances in accordance with Rumbell/Traub models for the axon.
                 
Edit History: Created by Nilapratim Sengupta in September 2023.
"""
def proc_setChannelConductances_pyramidalAxon(pyramidalCell, type_cell, index_parameter):
    
    # Dummy conductance parameters for ion channels
    dphi_cad   = 52 / 2e-3
    dbeta_cad  = 1 / 100	# In the paper beta = 50 [ms]

    dgnaf	= 0.2501575996 # Tuned NernstTest12 // 0.312 // tuned 84 // 150e-3 * 1.25 // default
    dgkdr  	= 0.0513842372 # Tuned NernstTest12 // 1.554 // tuned 84 // 125e-3
    dgka   	= 0
    dgk2   	= 0
    dgkm	= 0.0016943554 # Tuned NernstTest12 // 0.00138 // tuned 84 // 2.5 * 1.5e-3 * 2
    dgkc	= 0
    dgkahp	= 0

    dgcad 	= 0
    dgnap 	= 0
    dgcat 	= 0
    dgcal 	= 0
    
    # Dummy parameters for passive channels
    if (type_cell == 2):
        dgar = 2.83162e-05 # 1e-10
        dg_pas = 3.95094e-05 # 9.0916e-05
        de_pas = -67.92914 # -69.9598
        dcm = 1.09243 # 1.0 // (manual) // 0.9 // (default)
    elif (type_cell == 3):
        dgar = 3.61324e-05 # 5.0560e-05
        dg_pas = 1.71951e-05 # 1.7479e-05
        de_pas = -76.7721 # -79.9509
        dcm = 0.70449 # 1.0 // (manual) // 0.9 // (default)

    dRa	  	= 150	# (manual) // 250 // (default)
    dena	= 65.204611 # Using nernst( 115 , 5 , 1 )  // 50 // (default)
    dek		= -75.497 # Using nernst( 115 , 5 , 1 ) // -95 // (default)
    
    # Acquiring the list of ion channel parameter values based on type of cell and parameter index
    list_channelParameters = func_loadChannelParameters(type_cell, index_parameter)
    
    # Assigning values to designated parameters
    dgnaf = list_channelParameters[0]
    dgkdr = list_channelParameters[1]
    dgkm = list_channelParameters[2]
    dgka = list_channelParameters[3]
    dgkahp = list_channelParameters[4]
    dgnap = list_channelParameters[5]
    dgcal = list_channelParameters[6]
    dphi_cad = list_channelParameters[7]
    dbeta_cad = list_channelParameters[8]
    dtmmnaf = list_channelParameters[9]
    dtmhnaf = list_channelParameters[10]
    dtmkdr = list_channelParameters[11]
    dtmkm = list_channelParameters[12]
    dtmmka = list_channelParameters[13]
    dtmhka = list_channelParameters[14]
    dtmnap = list_channelParameters[15]
    dtmcal = list_channelParameters[16]
    dvsnaf = list_channelParameters[17]
    dvskdr = list_channelParameters[18]
    dvskm = list_channelParameters[19]
    dvska = list_channelParameters[20]
    dvsnap = list_channelParameters[21]
    dvscal = list_channelParameters[22]

    # Assigning values to ion channel conductances for the axon
    for sec in pyramidalCell.pyramidalAxon.TotalAxon:
        sec.gbar_naf = dgnaf * 2.133
        sec.gbar_kdr = dgkdr * 3.2
        sec.gbar_ka = dgka / 15
        sec.gbar_k2 = dgk2
        sec.g_pas = dg_pas  # * 50, used in Traub
        sec.e_pas = de_pas
        sec.Ra = dRa  # * 0.4, used in Traub
        sec.cm = dcm
        sec.ena = dena
        sec.ek = dek
    
    for sec in pyramidalCell.pyramidalAxon.Nodes_axon:
        if any(mec in sec.psection()['density_mechs'] for mec in ['pas']):
            sec.g_pas = 0.02

    for sec in pyramidalCell.pyramidalAxon.MyelinatedAxon:
        if any(mec in sec.psection()['density_mechs'] for mec in ['naf']):
            sec.gbar_naf = 0
        if any(mec in sec.psection()['density_mechs'] for mec in ['kdr']):
            sec.gbar_kdr = 0
        if any(mec in sec.psection()['density_mechs'] for mec in ['ka']):
            sec.gbar_ka = 0
        if any(mec in sec.psection()['density_mechs'] for mec in ['k2']):
            sec.gbar_k2 = 0
        if any(mec in sec.psection()['density_mechs'] for mec in ['pas']):
            sec.g_pas = dg_pas
        sec.cm = dcm * 0.04
        
# End of proc_setChannelConductances_pyramidalAxon()





"""
Name: proc_setChannelConductances_pyramidalCollateral()
Type: Procedure

Input Parameter(s): 1. Object - pyramidalCell
                    2. Number of collaterals
                    3. Type of cell (1: default, 2: young; 3: aged)
                    4. Parameter index (Either 0, 1 or 2)
Return Parameter: None

Description: Procedure to set ion channel conductances in accordance with Rumbell/Traub models for the collateral(s).
                 
Edit History: Created by Nilapratim Sengupta in September 2023.
"""
def proc_setChannelConductances_pyramidalCollateral(pyramidalCell, count_collaterals, type_cell, index_parameter):
    
    # Dummy conductance parameters for ion channels
    dphi_cad   = 52 / 2e-3
    dbeta_cad  = 1 / 100	# In the paper beta = 50 [ms]

    dgnaf	= 0.2501575996 # Tuned NernstTest12 // 0.312 // tuned 84 // 150e-3 * 1.25 // default
    dgkdr  	= 0.0513842372 # Tuned NernstTest12 // 1.554 // tuned 84 // 125e-3
    dgka   	= 0
    dgk2   	= 0
    dgkm	= 0.0016943554 # Tuned NernstTest12 // 0.00138 // tuned 84 // 2.5 * 1.5e-3 * 2
    dgkc	= 0
    dgkahp	= 0

    dgcad 	= 0
    dgnap 	= 0
    dgcat 	= 0
    dgcal 	= 0
    
    # Dummy parameters for passive channels
    if (type_cell == 2):
        dgar = 2.83162e-05 # 1e-10
        dg_pas = 3.95094e-05 # 9.0916e-05
        de_pas = -67.92914 # -69.9598
        dcm = 1.09243 # 1.0 // (manual) // 0.9 // (default)
    elif (type_cell == 3):
        dgar = 3.61324e-05 # 5.0560e-05
        dg_pas = 1.71951e-05 # 1.7479e-05
        de_pas = -76.7721 # -79.9509
        dcm = 0.70449 # 1.0 // (manual) // 0.9 // (default)

    dRa	  	= 150	# (manual) // 250 // (default)
    dena	= 65.204611 # Using nernst( 115 , 5 , 1 )  // 50 // (default)
    dek		= -75.497 # Using nernst( 115 , 5 , 1 ) // -95 // (default)
    
    # Acquiring the list of ion channel parameter values based on type of cell and parameter index
    list_channelParameters = func_loadChannelParameters(type_cell, index_parameter)
    
    # Assigning values to designated parameters
    dgnaf = list_channelParameters[0]
    dgkdr = list_channelParameters[1]
    dgkm = list_channelParameters[2]
    dgka = list_channelParameters[3]
    dgkahp = list_channelParameters[4]
    dgnap = list_channelParameters[5]
    dgcal = list_channelParameters[6]
    dphi_cad = list_channelParameters[7]
    dbeta_cad = list_channelParameters[8]
    dtmmnaf = list_channelParameters[9]
    dtmhnaf = list_channelParameters[10]
    dtmkdr = list_channelParameters[11]
    dtmkm = list_channelParameters[12]
    dtmmka = list_channelParameters[13]
    dtmhka = list_channelParameters[14]
    dtmnap = list_channelParameters[15]
    dtmcal = list_channelParameters[16]
    dvsnaf = list_channelParameters[17]
    dvskdr = list_channelParameters[18]
    dvskm = list_channelParameters[19]
    dvska = list_channelParameters[20]
    dvsnap = list_channelParameters[21]
    dvscal = list_channelParameters[22]
        
    # Assigning values to ion channel conductances for the collateral(s)
    for loopCounter in range(count_collaterals):
        for sec in pyramidalCell.pyramidalCollateral[loopCounter].TotalCollateral:
            sec.gbar_naf = dgnaf * 2.133
            sec.gbar_kdr = dgkdr * 3.2
            sec.gbar_ka = dgka / 15
            sec.gbar_k2 = dgk2
            sec.g_pas = dg_pas  # * 50, used in Traub
            sec.e_pas = de_pas
            sec.Ra = dRa  # * 0.4, used in Traub
            sec.cm = dcm
            sec.ena = dena
            sec.ek = dek

        for sec in pyramidalCell.pyramidalCollateral[loopCounter].Nodes_collateral:
            if any(mec in sec.psection()['density_mechs'] for mec in ['pas']):
                sec.g_pas = 0.02

        for sec in pyramidalCell.pyramidalCollateral[loopCounter].MyelinatedCollateral:
            if any(mec in sec.psection()['density_mechs'] for mec in ['naf']):
                sec.gbar_naf = 0
            if any(mec in sec.psection()['density_mechs'] for mec in ['kdr']):
                sec.gbar_kdr = 0
            if any(mec in sec.psection()['density_mechs'] for mec in ['ka']):
                sec.gbar_ka = 0
            if any(mec in sec.psection()['density_mechs'] for mec in ['k2']):
                sec.gbar_k2 = 0
            if any(mec in sec.psection()['density_mechs'] for mec in ['pas']):
                sec.g_pas = dg_pas
            sec.cm = dcm * 0.04           

# End of proc_setChannelConductances_pyramidalCollateral()





"""
Name: proc_setChannelKinetics()
Type: Procedure

Input Parameter(s): 1. Object - pyramidalCell
                    2. Type of cell (1: default, 2: young; 3: aged)
                    3. Parameter index (Either 0, 1 or 2)
Return Parameter: None

Description: Procedure to set ion channel kinetics in accordance with Rumbell/Traub models.
                 
Edit History: Created by Nilapratim Sengupta in July-August 2023.
"""
def proc_setChannelKinetics(pyramidalCell, type_cell, index_parameter):

    # Dummy kinetics parameters (taumods) for ion channels
    dtmmnaf = 1  # default
    dtmhnaf = 1  # default
    dtmkdr = 1   # default
    dtmkm = 1    # default
    dtmkc = 1    # default
    dtmcal = 1   # default
    dtmmka = 1   # default
    dtmhka = 1   # default
    dtmnap = 1   # default
    dtmar = 1    # default

    # Dummy kinetics parameters (vshifts) for ion channels
    dvsnaf = -3.5   # default
    dvskdr = 0      # default
    dvskm = 0       # default
    dvskc = 0       # default
    dvscal = 0      # default
    dvska = 0       # default
    dvsnap = 0      # default
    dvsar = 0       # default
    
    # Acquiring the list of ion channel parameter values based on type of cell and parameter index
    list_channelParameters = func_loadChannelParameters(type_cell, index_parameter)
    
    # Assigning values to designated parameters
    dgnaf = list_channelParameters[0]
    dgkdr = list_channelParameters[1]
    dgkm = list_channelParameters[2]
    dgka = list_channelParameters[3]
    dgkahp = list_channelParameters[4]
    dgnap = list_channelParameters[5]
    dgcal = list_channelParameters[6]
    dphi_cad = list_channelParameters[7]
    dbeta_cad = list_channelParameters[8]
    dtmmnaf = list_channelParameters[9]
    dtmhnaf = list_channelParameters[10]
    dtmkdr = list_channelParameters[11]
    dtmkm = list_channelParameters[12]
    dtmmka = list_channelParameters[13]
    dtmhka = list_channelParameters[14]
    dtmnap = list_channelParameters[15]
    dtmcal = list_channelParameters[16]
    dvsnaf = list_channelParameters[17]
    dvskdr = list_channelParameters[18]
    dvskm = list_channelParameters[19]
    dvska = list_channelParameters[20]
    dvsnap = list_channelParameters[21]
    dvscal = list_channelParameters[22]

    # Setting taumod values to the optimiser values
    h.mtaumod_naf = dtmmnaf
    h.htaumod_naf = dtmhnaf
    h.taumod_kdr = dtmkdr
    h.taumod_km = dtmkm
    h.taumod_kc = dtmkc
    h.taumod_cal = dtmcal
    h.mtaumod_ka = dtmmka
    h.htaumod_ka = dtmhka
    h.taumod_nap = dtmnap
    h.taumod_ar = dtmar

    # Setting usetable values to 0 if needed
    if h.mtaumod_naf != 1:
        h.usetable_naf = 0
    if h.htaumod_naf != 1:
        h.usetable_naf = 0
    if h.taumod_kdr != 1:
        h.usetable_kdr = 0
    if h.taumod_km != 1:
        h.usetable_km = 0
    if h.taumod_kc != 1:
        h.usetable_kc = 0
    if h.taumod_cal != 1:
        h.usetable_cal = 0
    if h.mtaumod_ka != 1:
        h.usetable_ka = 0
    if h.htaumod_ka != 1:
        h.usetable_ka = 0
    if h.taumod_nap != 1:
        h.usetable_nap = 0
    if h.taumod_ar != 1:
        h.usetable_ar = 0

    # Setting values for vshifts
    h.fastNashift_naf = dvsnaf
    h.vshift_kdr = dvskdr
    h.vshift_km = dvskm
    h.vshift_kc = dvskc
    h.vshift_cal = dvscal
    h.vshift_ka = dvska
    h.vshift_nap = dvsnap
    h.vshift_ar = dvsar

    # Setting usetable values to 0 if needed
    if h.fastNashift_naf != -3.5:
        h.usetable_naf = 0
    if h.vshift_kdr != 0:
       h.usetable_kdr = 0
    if h.vshift_km != 0:
        h.usetable_km = 0
    if h.vshift_kc != 0:
        h.usetable_kc = 0
    if h.vshift_cal != 0:
        h.usetable_cal = 0
    if h.vshift_ka != 0:
        h.usetable_ka = 0
    if h.vshift_nap != 0:
        h.usetable_nap = 0
    if h.vshift_ar != 0:
        h.usetable_ar = 0

# End of proc_setChannelKinetics()





"""
Name: proc_setParameters_final()
Type: Procedure

Input Parameter(s): 1. Object - pyramidalCell
                    2. Number of collaterals
                    3. Scale factor for gPas
                    4. Scale factor for gNa
                    5. Scale factor for gK
Return Parameter: None

Description: Procedure to update parameters just before run().
             This procedure helps adjust parameters for tuning.
             The hard-coded scale factors were determined from I-V curve studies while tuning the ion channels in accordance with the Scurfield-Latimer model.

                 
Edit History: Created by Nilapratim Sengupta in August 2023.
"""
def proc_setParameters_final(pyramidalCell, count_collaterals, scaleFactor_gPas, scaleFactor_gNa, scaleFactor_gK):
    
    # For the axon
    for sec in pyramidalCell.pyramidalAxon.Nodes_axon:
        sec.Ra = 150
        sec.cm = 1
        sec.g_pas = scaleFactor_gPas * 0.156000
        sec.gbar_naf = scaleFactor_gNa * 6
        sec.gbar_kdr = scaleFactor_gK * 0.0279
        # sec.gbar_ka = 8.33e-05
        # sec.gbar_k2 = 0

    for sec in pyramidalCell.pyramidalAxon.Paranodes_axon:
        sec.Ra = 150
        sec.cm = 1
        sec.g_pas = scaleFactor_gPas * 0.0005
        sec.gbar_naf = scaleFactor_gNa * 0.006
        sec.gbar_kdr = scaleFactor_gK * 0
        # sec.gbar_ka = 0
        # sec.gbar_k2 = 0

    for sec in pyramidalCell.pyramidalAxon.Juxtaparanodes_axon:
        sec.Ra = 150
        sec.cm = 1
        sec.g_pas = scaleFactor_gPas * 0.005
        sec.gbar_naf = scaleFactor_gNa * 0.06
        sec.gbar_kdr = scaleFactor_gK * 0.279
        # sec.gbar_ka = 0
        # sec.gbar_k2 = 0

    for sec in pyramidalCell.pyramidalAxon.Internodes_axon:
        sec.Ra = 150
        sec.cm = 1
        sec.g_pas = scaleFactor_gPas * 0.005
        sec.gbar_naf = scaleFactor_gNa * 0.06
        sec.gbar_kdr = scaleFactor_gK * 0.093
        # sec.gbar_ka = 0
        # sec.gbar_k2 = 0

    # For the collateral(s)
    if (count_collaterals != 0):
        for loopCounter in range(count_collaterals):
            for sec in pyramidalCell.pyramidalCollateral[loopCounter].Nodes_collateral:
                sec.Ra = 150
                sec.cm = 1
                sec.g_pas = scaleFactor_gPas * 0.156000
                sec.gbar_naf = scaleFactor_gNa * 6
                sec.gbar_kdr = scaleFactor_gK * 0.0279
                # sec.gbar_ka = 8.33e-05
                # sec.gbar_k2 = 0

            for sec in pyramidalCell.pyramidalCollateral[loopCounter].Paranodes_collateral:
                sec.Ra = 150
                sec.cm = 1
                sec.g_pas = scaleFactor_gPas * 0.0005
                sec.gbar_naf = scaleFactor_gNa * 0.006
                sec.gbar_kdr = scaleFactor_gK * 0
                # sec.gbar_ka = 0
                # sec.gbar_k2 = 0

            for sec in pyramidalCell.pyramidalCollateral[loopCounter].Juxtaparanodes_collateral:
                sec.Ra = 150
                sec.cm = 1
                sec.g_pas = scaleFactor_gPas * 0.005
                sec.gbar_naf = scaleFactor_gNa * 0.06
                sec.gbar_kdr = scaleFactor_gK * 0.279
                # sec.gbar_ka = 0
                # sec.gbar_k2 = 0

            for sec in pyramidalCell.pyramidalCollateral[loopCounter].Internodes_collateral:
                sec.Ra = 150
                sec.cm = 1
                sec.g_pas = scaleFactor_gPas * 0.005
                sec.gbar_naf = scaleFactor_gNa * 0.06
                sec.gbar_kdr = scaleFactor_gK * 0.093
                # sec.gbar_ka = 0
                # sec.gbar_k2 = 0
            
# End of proc_setParameters_final()





"""
Name: func_setRecLocations_pyramidalSomaDendrites()
Type: Function

Input Parameter(s): 1. Object - pyramidalCell
Return Parameter: List of recording locations (segments) 

Description: Function to set up recording locations in the soma and dendrite(s).

Edit History: Created by Nilapratim Sengupta in September 2023.
"""
def func_setRecLocations_pyramidalSomaDendrites(pyramidalCell):
    
    # Creating list for storing recording locations in the soma and dendrites
    list_recLocations_somaDendrites = []
    
    # Appending recording locations to the list
    list_recLocations_somaDendrites.append(pyramidalCell.pyramidalSomaDendrites.comp[1](0.5)) # soma
    list_recLocations_somaDendrites.append(pyramidalCell.pyramidalSomaDendrites.comp[68](0.5)) # distal apical dendrite
    
    # Returning the list of vectors 
    return list_recLocations_somaDendrites

# End of func_setRecLocations_pyramidalSomaDendrites()





"""
Name: func_setRecLocations_pyramidalAxon()
Type: Function

Input Parameter(s): 1. Object - pyramidalCell
                    2. Number of myelinated segments in the axon
Return Parameter: List of recording locations (segments) 

Description: Function to set up recording locations along the axon.

Edit History: Created by Nilapratim Sengupta in September 2023.
"""
def func_setRecLocations_pyramidalAxon(pyramidalCell, count_axonSegments):
    
    # Creating list for storing recording locations along the axon
    list_recLocations_axon = []
    
    # Recording locations - for the axon
    location0_axon = 0
    location1_axon = math.floor(count_axonSegments / 2)
    location2_axon = (count_axonSegments - 1) # Penultimate: because of current reflection in the ultimate node
    
    # Appending recording locations to the list
    list_recLocations_axon.append(pyramidalCell.pyramidalAxon.hill_axon(0.5))
    list_recLocations_axon.append(pyramidalCell.pyramidalAxon.iseg_axon(0.5))
    list_recLocations_axon.append(pyramidalCell.pyramidalAxon.nodes_axon[location0_axon](0.5))
    list_recLocations_axon.append(pyramidalCell.pyramidalAxon.nodes_axon[location1_axon](0.5))
    list_recLocations_axon.append(pyramidalCell.pyramidalAxon.nodes_axon[location2_axon](0.5))
    list_recLocations_axon.append(pyramidalCell.pyramidalAxon.paranodeThrees_axon[2*location0_axon](0.5))
    list_recLocations_axon.append(pyramidalCell.pyramidalAxon.paranodeThrees_axon[2*location1_axon](0.5))
    list_recLocations_axon.append(pyramidalCell.pyramidalAxon.paranodeThrees_axon[2*location2_axon](0.5))
    list_recLocations_axon.append(pyramidalCell.pyramidalAxon.juxtaparanodes_axon[2*location0_axon](0.5))
    list_recLocations_axon.append(pyramidalCell.pyramidalAxon.juxtaparanodes_axon[2*location1_axon](0.5))
    list_recLocations_axon.append(pyramidalCell.pyramidalAxon.juxtaparanodes_axon[2*location2_axon](0.5))
    list_recLocations_axon.append(pyramidalCell.pyramidalAxon.internodes_axon[location0_axon](0.5))
    list_recLocations_axon.append(pyramidalCell.pyramidalAxon.internodes_axon[location1_axon](0.5))
    list_recLocations_axon.append(pyramidalCell.pyramidalAxon.internodes_axon[location2_axon](0.5))
    
    # Returning the list of vectors 
    return list_recLocations_axon

# End of func_setRecLocations_pyramidalAxon()





"""
Name: func_setRecLocations_pyramidalCollateral()
Type: Function

Input Parameter(s): 1. Object - pyramidalCell
                    2. Number of collaterals
                    3. List of number of myelinated segments in the collateral(s)
Return Parameter: List of recording locations (segments) 

Description: Function to set up recording locations along the collateral(s).

Edit History: Created by Nilapratim Sengupta in September 2023.
"""
def func_setRecLocations_pyramidalCollateral(pyramidalCell, count_collaterals, countList_collateralSegments):
    
    # Creating list for storing recording locations along the collateral(s)
    list_recLocations_collaterals = []
    
    # Looping over the collateral(s)
    for loopCounter_collateral in range(count_collaterals):
        
        # Recording locations - for the collateral(s)
        location0_collateral = 0
        location1_collateral = math.floor(countList_collateralSegments[loopCounter_collateral] / 2)
        location2_collateral = (countList_collateralSegments[loopCounter_collateral] - 1) # Penultimate: because of current reflection in the ultimate node 
        
        list_recLocations_collaterals.append(pyramidalCell.pyramidalCollateral[loopCounter_collateral].nodes_collateral[location0_collateral](0.5))
        list_recLocations_collaterals.append(pyramidalCell.pyramidalCollateral[loopCounter_collateral].nodes_collateral[location1_collateral](0.5))
        list_recLocations_collaterals.append(pyramidalCell.pyramidalCollateral[loopCounter_collateral].nodes_collateral[location2_collateral](0.5))
        list_recLocations_collaterals.append(pyramidalCell.pyramidalCollateral[loopCounter_collateral].paranodeThrees_collateral[2*location0_collateral](0.5))
        list_recLocations_collaterals.append(pyramidalCell.pyramidalCollateral[loopCounter_collateral].paranodeThrees_collateral[2*location1_collateral](0.5))
        list_recLocations_collaterals.append(pyramidalCell.pyramidalCollateral[loopCounter_collateral].paranodeThrees_collateral[2*location2_collateral](0.5))
        list_recLocations_collaterals.append(pyramidalCell.pyramidalCollateral[loopCounter_collateral].juxtaparanodes_collateral[2*location0_collateral](0.5))
        list_recLocations_collaterals.append(pyramidalCell.pyramidalCollateral[loopCounter_collateral].juxtaparanodes_collateral[2*location1_collateral](0.5))
        list_recLocations_collaterals.append(pyramidalCell.pyramidalCollateral[loopCounter_collateral].juxtaparanodes_collateral[2*location2_collateral](0.5))
        list_recLocations_collaterals.append(pyramidalCell.pyramidalCollateral[loopCounter_collateral].internodes_collateral[location0_collateral](0.5))
        list_recLocations_collaterals.append(pyramidalCell.pyramidalCollateral[loopCounter_collateral].internodes_collateral[location1_collateral](0.5))
        list_recLocations_collaterals.append(pyramidalCell.pyramidalCollateral[loopCounter_collateral].internodes_collateral[location2_collateral](0.5))
    
    # Returning the list of vectors 
    return list_recLocations_collaterals

# End of func_setRecLocations_pyramidalCollateral()
    




"""
Name: func_detectSpikeTimes_pyramidalSomaDendrites()
Type: Function

Input Parameter(s): 1. List of recording locations (segments) in the soma and dendrites
Return Parameter: List of vectors (HOC)

Description: Function to set up spike detectors in the soma and dendrite(s).

Edit History: Created by Nilapratim Sengupta in August-September 2023.
"""
def func_detectSpikeTimes_pyramidalSomaDendrites(list_recLocations_somaDendrites):
    
    # Creating netcon objects at the recording locations  
    netCons_somaDendrites = [h.NetCon(seg._ref_v, None, sec=seg.sec) for seg in list_recLocations_somaDendrites] # seg.sec = section in which the segment exists
    
    # Creating list of vectors to record 
    vecList_spikeTimes_somaDendrites = [h.Vector() for _ in netCons_somaDendrites]
    
    # Appending vectors to the list - note the order of vectors
    for vec_spikeTimes, netCons in zip(vecList_spikeTimes_somaDendrites, netCons_somaDendrites):
        netCons.record(vec_spikeTimes)
    
    # Returning the list of vectors 
    return vecList_spikeTimes_somaDendrites

# End of func_detectSpikeTimes_pyramidalSomaDendrites()





"""
Name: func_detectSpikeTimes_pyramidalAxon()
Type: Function

Input Parameter(s): 1. List of recording locations (segments) along the axon
Return Parameter: List of vectors (HOC) 

Description: Function to set up spike detectors at different locations along the axon.

Edit History: Created by Nilapratim Sengupta in August-September 2023.
"""
def func_detectSpikeTimes_pyramidalAxon(list_recLocations_axon):
    
    # Creating netcon objects at the recording locations  
    netCons_axon = [h.NetCon(seg._ref_v, None, sec=seg.sec) for seg in list_recLocations_axon] # seg.sec = section in which the segment exists
    
    # Creating list of vectors to record 
    vecList_spikeTimes_axon = [h.Vector() for _ in netCons_axon]
    
    # Appending vectors to the list - note the order of vectors
    for vec_spikeTimes, netCons in zip(vecList_spikeTimes_axon, netCons_axon):
        netCons.record(vec_spikeTimes)
    
    # Returning the list of vectors 
    return vecList_spikeTimes_axon

# End of func_detectSpikeTimes_pyramidalAxon()





"""
Name: func_detectSpikeTimes_pyramidalCollateral()
Type: Function

Input Parameter(s): 1. List of recording locations (segments) along the collateral(s)
Return Parameter: List of vectors (HOC) 

Description: Function to set up spike detectors at different locations along the collateral(s).

Edit History: Created by Nilapratim Sengupta in August-September 2023.
"""
def func_detectSpikeTimes_pyramidalCollateral(list_recLocations_collaterals):    
    
    # Creating netcon objects at the recording locations  
    netCons_collaterals = [h.NetCon(seg._ref_v, None, sec=seg.sec) for seg in list_recLocations_collaterals] # seg.sec = section in which the segment exists
    
    # Creating list of vectors to record 
    vecList_spikeTimes_collaterals = [h.Vector() for _ in netCons_collaterals] 
    
    # Appending vectors to the list - note the order of vectors
    for vec_spikeTimes, netCons in zip(vecList_spikeTimes_collaterals, netCons_collaterals):
        netCons.record(vec_spikeTimes)
    
    # Returning the list of vectors 
    return vecList_spikeTimes_collaterals

# End of func_detectSpikeTimes_pyramidalCollateral()





"""
Name: func_calculateFiringRate()
Type: Function

Input Parameter(s): 1. Vector with spike times
                    2. Start time
                    3. Stop time
Return Parameter: Firing rate 

Description: Function to calculate the firing rate within a time-window from a vector containing spike times.

Edit History: Created by Nilapratim Sengupta in August 2023.
"""
def func_calculateFiringRate(vec_spikeTimes, time_start, time_stop):
    
    # Creating lists to store ISI and FR
    list_interSpikeInterval = []
    list_firingRate = []
    
    # Calculating FR from ISI and appending to lists
    for loopCounter in range(len(vec_spikeTimes) - 1):
        if vec_spikeTimes[loopCounter] >= time_start and vec_spikeTimes[loopCounter + 1] <= time_stop:
            interSpikeInterval = vec_spikeTimes[loopCounter + 1] - vec_spikeTimes[loopCounter]
            list_interSpikeInterval.append(interSpikeInterval)
            firingRate = 1000 / interSpikeInterval
            list_firingRate.append(firingRate)

    # Returning FR
    if len(list_firingRate) == 0:
        return 0.0
    else:
        return np.mean(list_firingRate)
    
# End of func_calculateFiringRate()





"""
Name: func_calculateAPCount()
Type: Function

Input Parameter(s): 1. Vector with spike times
                    2. Start time
                    3. Stop time
Return Parameter: AP count 

Description: Function to calculate the AP count within a time-window from a vector containing spike times.

Edit History: Created by Nilapratim Sengupta in August 2023.
"""
def func_calculateAPCount(vec_spikeTimes, time_start, time_stop):
    
    # Calculating AP count
    countAP = 0
    for loopCounter in range(len(vec_spikeTimes)):
        if vec_spikeTimes[loopCounter] >= time_start and vec_spikeTimes[loopCounter] <= time_stop:
            countAP = (countAP + 1)

    # Returning AP count
    return countAP
    
# End of func_calculateAPCount()





"""
Name: func_calculateConductionTime()
Type: Function

Input Parameter(s): 1. Vector with spike times at start location
                    2. Vector with spike times at stop location
                    2. Start time
                    3. Stop time
Return Parameter: Conduction time

Description: Function to calculate the firing rate and AP count within a time-window from a vector containing spike times.

Edit History: Created by Nilapratim Sengupta in August 2023.
              Edited by Nilapratim Sengupta in May 2024 to focus on corresponding spikes.
"""
def func_calculateConductionTime(spikeTimes_startLocation, spikeTimes_stopLocation, time_start, time_stop):
    
    # Creating list to store conduction times
    list_conductionTime = []
    
    # Calculating APCs at the start and stop locations
    countAP_startLocation = func_calculateAPCount(spikeTimes_startLocation, time_start, time_stop)
    countAP_stopLocation = func_calculateAPCount(spikeTimes_stopLocation, time_start, time_stop)
    
    # Calculating conduction time
    if countAP_startLocation > 0 and countAP_stopLocation > 0:
        for loopCounter in range(countAP_stopLocation):
            list_conductionTime.append(spikeTimes_stopLocation[loopCounter] - spikeTimes_startLocation[loopCounter])
            # Displaying on the console
            print(f'\nConduction Time for spike {loopCounter+1} = {spikeTimes_stopLocation[loopCounter] - spikeTimes_startLocation[loopCounter]}') # Console-display
        conductionTime = np.mean(list_conductionTime)
    else:
        conductionTime = 0

    # Returning conduction time
    return conductionTime
  
# End of func_calculateConductionTime()





"""
Name: func_calculatePropagationDetails_pyramidalSomaDendrites()
Type: Function

Input Parameter(s): 1. Object - pyramidalCell
                    2. List of vectors with spike times from the soma and dendrites
                    3. Start time
                    4. Stop time
Return Parameter: List - 1. Firing rate at start location
                         2. Firing rate at stop location
                         3. AP count at start location
                         4. AP count at stop location
                         5. Percentage of APs lost (between start and stop locations)
                         6. Conduction velocity (between start and stop locations)

Description: Function to calculate details of AP propagation in the soma and dendrites.

Edit History: Created by Nilapratim Sengupta in August 2023.
"""
def func_calculatePropagationDetails_pyramidalSomaDendrites(pyramidalCell, vecList_spikeTimes_somaDendrites, time_start, time_stop):

    # Calculating FRs at the start and stop locations
    firingRate_startLocation = func_calculateFiringRate(vecList_spikeTimes_somaDendrites[0], time_start, time_stop)
    firingRate_stopLocation = func_calculateFiringRate(vecList_spikeTimes_somaDendrites[1], time_start, time_stop)
    
    # Calculating APCs at the start and stop locations
    countAP_startLocation = func_calculateAPCount(vecList_spikeTimes_somaDendrites[0], time_start, time_stop)
    countAP_stopLocation = func_calculateAPCount(vecList_spikeTimes_somaDendrites[1], time_start, time_stop)
    
    # Calculating percentage of APs lost during propagation (between start and stop locations)
    if (countAP_startLocation != 0):
        percentage_lostAP = ((countAP_startLocation - countAP_stopLocation) / countAP_startLocation) * 100
        if (percentage_lostAP < 0):
            percentage_lostAP = 0 # When more spikes are detected at stop location
    else: 
        percentage_lostAP = 100 # When no spikes detected at start location
        
    # Calculating conduction velocity between start and stop locations
    conductionTime = func_calculateConductionTime(vecList_spikeTimes_somaDendrites[0], vecList_spikeTimes_somaDendrites[1], time_start, time_stop)
    conductionDistance = h.distance(pyramidalCell.pyramidalSomaDendrites.comp[1](0.5), pyramidalCell.pyramidalSomaDendrites.comp[68](0.5))
    if conductionTime != 0:
        conductionVelocity = (1e-3 * (conductionDistance / conductionTime))
    else:
        conductionVelocity = 0
    
    # Determining outcome flag for propagation
    outcomeFlag = 0
    if (countAP_startLocation != 0 and countAP_stopLocation != 0):
        outcomeFlag = 1
    
    # Returning list of parameters
    list_propagationDetails_somaDendrites = []
    list_propagationDetails_somaDendrites.append(firingRate_startLocation)
    list_propagationDetails_somaDendrites.append(firingRate_stopLocation)
    list_propagationDetails_somaDendrites.append(countAP_startLocation)
    list_propagationDetails_somaDendrites.append(countAP_stopLocation)
    list_propagationDetails_somaDendrites.append(percentage_lostAP)
    list_propagationDetails_somaDendrites.append(conductionVelocity)
    list_propagationDetails_somaDendrites.append(outcomeFlag)
    
    return list_propagationDetails_somaDendrites
    
# End of func_calculatePropagationDetails_pyramidalSomaDendrites()





"""
Name: func_calculatePropagationDetails_pyramidalAxon()
Type: Function

Input Parameter(s): 1. Object - pyramidalCell
                    2. List of vectors with spike times along the axon
                    3. Number of myelinated segments in the axon
                    4. Start location (0 or 1)
                    5. Stop location (1 or 2)
                    6. Start time
                    7. Stop time
Return Parameter: List - 1. Firing rate at start location
                         2. Firing rate at stop location
                         3. AP count at start location
                         4. AP count at stop location
                         5. Percentage of APs lost (between start and stop locations)
                         6. Conduction velocity (between start and stop locations)
                         7. Outcome flag (depends on signal being active/passive in designated regions) - 0: Not Acceptable, 1: Acceptable

Description: Function to calculate details of AP propagation along the axon.

Edit History: Created by Nilapratim Sengupta in August 2023.
"""
def func_calculatePropagationDetails_pyramidalAxon(pyramidalCell, vecList_spikeTimes_axon, count_axonSegments, location_start_axon, location_stop_axon, time_start, time_stop):

    # Recording locations - for the axon
    location0_axon = 0
    location1_axon = math.floor(count_axonSegments / 2)
    location2_axon = (count_axonSegments - 1) # Penultimate: because of current reflection in the ultimate node
    
    # Obtaining spikeTime vectors for hillock and iseg
    spikeTimes_startLocation_hill_axon = vecList_spikeTimes_axon[0]
    spikeTimes_startLocation_iseg_axon = vecList_spikeTimes_axon[1]
    
    # Selecting spikeTime vectors for start location
    if (location_start_axon == 0):
        startLocation_node_axon = pyramidalCell.pyramidalAxon.nodes_axon[location0_axon](0.5)
        spikeTimes_startLocation_node_axon = vecList_spikeTimes_axon[2]
        spikeTimes_startLocation_paranode_axon = vecList_spikeTimes_axon[5]
        spikeTimes_startLocation_juxtaparanode_axon = vecList_spikeTimes_axon[8]
        spikeTimes_startLocation_internode_axon = vecList_spikeTimes_axon[11]
    elif (location_start_axon == 1):
        startLocation_node_axon = pyramidalCell.pyramidalAxon.nodes_axon[location1_axon](0.5)
        spikeTimes_startLocation_node_axon = vecList_spikeTimes_axon[3]
        spikeTimes_startLocation_paranode_axon = vecList_spikeTimes_axon[6]
        spikeTimes_startLocation_juxtaparanode_axon = vecList_spikeTimes_axon[9]
        spikeTimes_startLocation_internode_axon = vecList_spikeTimes_axon[12]
        
    # Selecting spikeTime vectors for stop location
    if (location_stop_axon == 1):
        stopLocation_node_axon = pyramidalCell.pyramidalAxon.nodes_axon[location1_axon](0.5)
        spikeTimes_stopLocation_node_axon = vecList_spikeTimes_axon[3]
        spikeTimes_stopLocation_paranode_axon = vecList_spikeTimes_axon[6]
        spikeTimes_stopLocation_juxtaparanode_axon = vecList_spikeTimes_axon[9]
        spikeTimes_stopLocation_internode_axon = vecList_spikeTimes_axon[12]
    elif (location_stop_axon == 2):
        stopLocation_node_axon = pyramidalCell.pyramidalAxon.nodes_axon[location2_axon](0.5)
        spikeTimes_stopLocation_node_axon = vecList_spikeTimes_axon[4]
        spikeTimes_stopLocation_paranode_axon = vecList_spikeTimes_axon[7]
        spikeTimes_stopLocation_juxtaparanode_axon = vecList_spikeTimes_axon[10]
        spikeTimes_stopLocation_internode_axon = vecList_spikeTimes_axon[13]

    # Calculating FRs at the start and stop locations
    firingRate_startLocation = func_calculateFiringRate(spikeTimes_startLocation_node_axon, time_start, time_stop)
    firingRate_stopLocation = func_calculateFiringRate(spikeTimes_stopLocation_node_axon, time_start, time_stop)
    
    # Calculating APCs at the start and stop locations
    countAP_startLocation_node = func_calculateAPCount(spikeTimes_startLocation_node_axon, time_start, time_stop)
    countAP_stopLocation_node = func_calculateAPCount(spikeTimes_stopLocation_node_axon, time_start, time_stop)
    
    # Calculating percentage of APs lost during propagation (between start and stop locations)
    if (countAP_startLocation_node != 0):
        percentage_lostAP = ((countAP_startLocation_node - countAP_stopLocation_node) / countAP_startLocation_node) * 100
        if (percentage_lostAP < 0):
            percentage_lostAP = 0 # When more spikes are detected at stop location
    else: 
        percentage_lostAP = 100 # When no spikes detected at start location
        
    # Calculating conduction velocity between start and stop locations
    conductionTime = func_calculateConductionTime(spikeTimes_startLocation_node_axon, spikeTimes_stopLocation_node_axon, time_start, time_stop)
    conductionDistance = h.distance(startLocation_node_axon, stopLocation_node_axon)
    # Display for the user
    print(f'\nFor Axon:\tConduction Distance: {conductionDistance}\tConduction Time: {conductionTime}') # Display for the user
    if conductionTime != 0:
        conductionVelocity = (1e-3 * (conductionDistance / conductionTime))
    else:
        conductionVelocity = 0
    
    # Determining outcome flag for propagation
    outcomeFlag = 0
    countAP_startLocation_paranode = func_calculateAPCount(spikeTimes_startLocation_paranode_axon, time_start, time_stop)
    countAP_startLocation_juxtaparanode = func_calculateAPCount(spikeTimes_startLocation_juxtaparanode_axon, time_start, time_stop)
    countAP_startLocation_internode = func_calculateAPCount(spikeTimes_startLocation_internode_axon, time_start, time_stop)
    countAP_stopLocation_paranode = func_calculateAPCount(spikeTimes_stopLocation_paranode_axon, time_start, time_stop)
    countAP_stopLocation_juxtaparanode = func_calculateAPCount(spikeTimes_stopLocation_juxtaparanode_axon, time_start, time_stop)
    countAP_stopLocation_internode = func_calculateAPCount(spikeTimes_stopLocation_internode_axon, time_start, time_stop)
    if (countAP_startLocation_node != 0 and 
        countAP_startLocation_paranode == 0 and
        countAP_startLocation_juxtaparanode == 0 and
        countAP_startLocation_internode == 0 and
        countAP_stopLocation_node != 0 and
        countAP_stopLocation_paranode == 0 and
        countAP_stopLocation_juxtaparanode == 0 and
        countAP_stopLocation_internode == 0):
        outcomeFlag = 1
    
    # Returning list of parameters
    list_propagationDetails_axon = []
    list_propagationDetails_axon.append(firingRate_startLocation)
    list_propagationDetails_axon.append(firingRate_stopLocation)
    list_propagationDetails_axon.append(countAP_startLocation_node)
    list_propagationDetails_axon.append(countAP_stopLocation_node)
    list_propagationDetails_axon.append(percentage_lostAP)
    list_propagationDetails_axon.append(conductionVelocity)
    list_propagationDetails_axon.append(outcomeFlag)
    
    return list_propagationDetails_axon
    
# End of func_calculatePropagationDetails_pyramidalAxon()





"""
Name: func_calculatePropagationDetails_pyramidalCollateral()
Type: Function

Input Parameter(s): 1. Object - pyramidalCell
                    2. List of vectors with spike times along the axon
                    3. List of vectors with spike times along the collateral(s)
                    4. Index of particular collateral
                    5. Number of myelinated segments in the particular collateral
                    6. Start location (-1 or 0 or 1)
                    7. Stop location (1 or 2)
                    8. Start time
                    9. Stop time
Return Parameter: List - 1. Firing rate at start location
                         2. Firing rate at stop location
                         3. AP count at start location
                         4. AP count at stop location
                         5. Percentage of APs lost (between start and stop locations)
                         6. Conduction velocity (between start and stop locations)
                         7. Outcome flag (depends on signal being active/passive in designated regions) - 0: Not Acceptable, 1: Acceptable

Description: Function to calculate details of AP propagation along the particular collateral.

Edit History: Created by Nilapratim Sengupta in August 2023.
              Edited by Nilapratim Sengupta in December 2023 to update options for location_start_collateral.
"""
def func_calculatePropagationDetails_pyramidalCollateral(pyramidalCell, vecList_spikeTimes_axon, vecList_spikeTimes_collaterals, index_collateral, count_collateralSegments, location_start_collateral, location_stop_collateral, time_start, time_stop):

    # Recording locations - for the collateral
    location0_axon = 0
    location0_collateral = 0
    location1_collateral = math.floor(count_collateralSegments / 2)
    location2_collateral = (count_collateralSegments - 1) # Penultimate: because of current reflection in the ultimate node
    
    # Selecting spikeTime vectors for start location - IMPORTANT: 12 denotes the number of recording locations along each collateral
    if (location_start_collateral == -1):
        startLocation_node_collateral = pyramidalCell.pyramidalAxon.nodes_axon[location0_axon](0.5)
        spikeTimes_startLocation_node_collateral = vecList_spikeTimes_axon[2]
        spikeTimes_startLocation_paranode_collateral = vecList_spikeTimes_axon[5]
        spikeTimes_startLocation_juxtaparanode_collateral = vecList_spikeTimes_axon[8]
        spikeTimes_startLocation_internode_collateral = vecList_spikeTimes_axon[11]
    elif (location_start_collateral == 0):
        startLocation_node_collateral = pyramidalCell.pyramidalCollateral[index_collateral].nodes_collateral[location0_collateral](0.5)
        spikeTimes_startLocation_node_collateral = vecList_spikeTimes_collaterals[12*index_collateral+0]
        spikeTimes_startLocation_paranode_collateral = vecList_spikeTimes_collaterals[12*index_collateral+3]
        spikeTimes_startLocation_juxtaparanode_collateral = vecList_spikeTimes_collaterals[12*index_collateral+6]
        spikeTimes_startLocation_internode_collateral = vecList_spikeTimes_collaterals[12*index_collateral+9]
    elif (location_start_collateral == 1):
        startLocation_node_collateral = pyramidalCell.pyramidalCollateral[index_collateral].nodes_collateral[location1_collateral](0.5)
        spikeTimes_startLocation_node_collateral = vecList_spikeTimes_collaterals[12*index_collateral+1]
        spikeTimes_startLocation_paranode_collateral = vecList_spikeTimes_collaterals[12*index_collateral+4]
        spikeTimes_startLocation_juxtaparanode_collateral = vecList_spikeTimes_collaterals[12*index_collateral+7]
        spikeTimes_startLocation_internode_collateral = vecList_spikeTimes_collaterals[12*index_collateral+10]
        
    # Selecting spikeTime vectors for stop location
    if (location_stop_collateral == 1):
        stopLocation_node_collateral = pyramidalCell.pyramidalCollateral[index_collateral].nodes_collateral[location1_collateral](0.5)
        spikeTimes_stopLocation_node_collateral = vecList_spikeTimes_collaterals[12*index_collateral+1]
        spikeTimes_stopLocation_paranode_collateral = vecList_spikeTimes_collaterals[12*index_collateral+4]
        spikeTimes_stopLocation_juxtaparanode_collateral = vecList_spikeTimes_collaterals[12*index_collateral+7]
        spikeTimes_stopLocation_internode_collateral = vecList_spikeTimes_collaterals[12*index_collateral+10]
    elif (location_stop_collateral == 2):
        stopLocation_node_collateral = pyramidalCell.pyramidalCollateral[index_collateral].nodes_collateral[location2_collateral](0.5)
        spikeTimes_stopLocation_node_collateral = vecList_spikeTimes_collaterals[12*index_collateral+2]
        spikeTimes_stopLocation_paranode_collateral = vecList_spikeTimes_collaterals[12*index_collateral+5]
        spikeTimes_stopLocation_juxtaparanode_collateral = vecList_spikeTimes_collaterals[12*index_collateral+8]
        spikeTimes_stopLocation_internode_collateral = vecList_spikeTimes_collaterals[12*index_collateral+11]

    # Calculating FRs at the start and stop locations
    firingRate_startLocation = func_calculateFiringRate(spikeTimes_startLocation_node_collateral, time_start, time_stop)
    firingRate_stopLocation = func_calculateFiringRate(spikeTimes_stopLocation_node_collateral, time_start, time_stop)
    
    # Calculating APCs at the start and stop locations
    countAP_startLocation_node = func_calculateAPCount(spikeTimes_startLocation_node_collateral, time_start, time_stop)
    countAP_stopLocation_node = func_calculateAPCount(spikeTimes_stopLocation_node_collateral, time_start, time_stop)
    
    # Calculating percentage of APs lost during propagation (between start and stop locations)
    if (countAP_startLocation_node != 0):
        percentage_lostAP = ((countAP_startLocation_node - countAP_stopLocation_node) / countAP_startLocation_node) * 100
        if (percentage_lostAP < 0):
            percentage_lostAP = 0 # When more spikes are detected at stop location
    else: 
        percentage_lostAP = 100 # When no spikes detected at start location
        
    # Calculating conduction velocity between start and stop locations
    conductionTime = func_calculateConductionTime(spikeTimes_startLocation_node_collateral, spikeTimes_stopLocation_node_collateral, time_start, time_stop)
    conductionDistance = h.distance(startLocation_node_collateral, stopLocation_node_collateral)
    # Display for the user
    print(f'\nFor Collateral: {index_collateral}\tConduction Distance: {conductionDistance}\tConduction Time: {conductionTime}') # Display for the user
    if conductionTime != 0:
        conductionVelocity = (1e-3 * (conductionDistance / conductionTime))
    else:
        conductionVelocity = 0
    
    # Determining outcome flag for propagation
    outcomeFlag = 0
    countAP_startLocation_paranode = func_calculateAPCount(spikeTimes_startLocation_paranode_collateral, time_start, time_stop)
    countAP_startLocation_juxtaparanode = func_calculateAPCount(spikeTimes_startLocation_juxtaparanode_collateral, time_start, time_stop)
    countAP_startLocation_internode = func_calculateAPCount(spikeTimes_startLocation_internode_collateral, time_start, time_stop)
    countAP_stopLocation_paranode = func_calculateAPCount(spikeTimes_stopLocation_paranode_collateral, time_start, time_stop)
    countAP_stopLocation_juxtaparanode = func_calculateAPCount(spikeTimes_stopLocation_juxtaparanode_collateral, time_start, time_stop)
    countAP_stopLocation_internode = func_calculateAPCount(spikeTimes_stopLocation_internode_collateral, time_start, time_stop)
    if (countAP_startLocation_node != 0 and 
        countAP_startLocation_paranode == 0 and
        countAP_startLocation_juxtaparanode == 0 and
        countAP_startLocation_internode == 0 and
        countAP_stopLocation_node != 0 and
        countAP_stopLocation_paranode == 0 and
        countAP_stopLocation_juxtaparanode == 0 and
        countAP_stopLocation_internode == 0):
        outcomeFlag = 1
    
    # Returning list of parameters
    list_propagationDetails_collateral = []
    list_propagationDetails_collateral.append(firingRate_startLocation)
    list_propagationDetails_collateral.append(firingRate_stopLocation)
    list_propagationDetails_collateral.append(countAP_startLocation_node)
    list_propagationDetails_collateral.append(countAP_stopLocation_node)
    list_propagationDetails_collateral.append(percentage_lostAP)
    list_propagationDetails_collateral.append(conductionVelocity)
    list_propagationDetails_collateral.append(outcomeFlag)
    
    return list_propagationDetails_collateral
    
# End of func_calculatePropagationDetails_pyramidalCollateral() 





"""
Name: func_recordVectors_pyramidalSomaDendrites()
Type: Function

Input Parameter(s): 1. List of recording locations (segments) in the soma and dendrites
Return Parameter: List of vectors (HOC)

Description: Function to set up vectors to record membrane potentials in the soma and dendrite(s).

Edit History: Created by Nilapratim Sengupta in August 2023.
"""
def func_recordVectors_pyramidalSomaDendrites(list_recLocations_somaDendrites):

    # Creating list of vectors to record 
    vecList_potentials_somaDendrites = [h.Vector().record(seg._ref_v) for seg in list_recLocations_somaDendrites] 
    
    # Returning list of vectors
    return vecList_potentials_somaDendrites

# End of func_recordVectors_pyramidalSomaDendrites() 





"""
Name: func_recordVectors_pyramidalAxon()
Type: Function

Input Parameter(s): 1. List of recording locations (segments) along the axon
Return Parameter: List of vectors (HOC) 

Description: Function to set up vectors to record membrane potentials along the axon.

Edit History: Created by Nilapratim Sengupta in August-September 2023.
"""
def func_recordVectors_pyramidalAxon(list_recLocations_axon):
    
    # Creating list of vectors to record 
    vecList_potentials_axon = [h.Vector().record(seg._ref_v) for seg in list_recLocations_axon] 
    
    # Returning the list of vectors 
    return vecList_potentials_axon

# End of func_recordVectors_pyramidalAxon()   





"""
Name: func_recordVectors_pyramidalCollateral()
Type: Function

Input Parameter(s): 1. List of recording locations (segments) along the collateral(s)
Return Parameter: List of vectors (HOC) 

Description: Function to set up vectors to record membrane potentials along the collateral(s).

Edit History: Created by Nilapratim Sengupta in September 2023.
"""
def func_recordVectors_pyramidalCollateral(list_recLocations_collaterals):
   
    # Creating list of vectors to record 
    vecList_potentials_collaterals = [h.Vector().record(seg._ref_v) for seg in list_recLocations_collaterals] 
    
    # Returning the list of vectors 
    return vecList_potentials_collaterals

# End of func_recordVectors_pyramidalCollateral()   





"""
Name: proc_plotPotentials_pyramidalSomaDendrites()
Type: Procedure

Input Parameter(s): 1. Object - matplotlib pyplot (plt)
                    2. List of recording locations (segments) in the soma and dendrites
                    3. List of membrane potential vectors for the soma and dendrites
                    4. Time vector
Return Parameter: None

Description: Procedure to plot membrane potential(s) in the soma and dendrite(s).

Edit History: Created by Nilapratim Sengupta in August-September 2023.
"""
def proc_plotPotentials_pyramidalSomaDendrites(plt, list_recLocations_somaDendrites, vecList_potentials_somaDendrites, vec_time):  

    # Creating trace and label lists
    list_traces_start = []
    list_traces_stop = []
    list_labels_start = []
    list_labels_stop = []
    
    # Appending traces to the lists of traces
    list_traces_start.append(vecList_potentials_somaDendrites[0])
    list_traces_stop.append(vecList_potentials_somaDendrites[1])
    
    # Appending corresponding labels to the lists of labels
    list_labels_start.append(list_recLocations_somaDendrites[0])
    list_labels_stop.append(list_recLocations_somaDendrites[1])

    # Creating a figure
    fig0, ax0 = plt.subplots(nrows=2, sharex=True, squeeze=False)
    ax0_1, ax0_2 = ax0.flatten()
    
    # Plotting in the figure
    for trace, label in zip(list_traces_start, list_labels_start):
        ax0_1.plot(vec_time, trace, label=label)
    for trace, label in zip(list_traces_stop, list_labels_stop):
        ax0_2.plot(vec_time, trace, label=label)

    # Adding figure details
    ax0_2.set_xlabel('Time (ms)')
    ax0_1.set_ylabel('Membrane Potential (mV)')
    ax0_2.set_ylabel('Membrane Potential (mV)')
    ax0_1.set_title('Soma')
    ax0_2.set_title('Dendrite')
    ax0_1.legend()
    ax0_2.legend()

    # Visualizing the figure
    # plt.show()
    
# End of proc_plotPotentials_pyramidalSomaDendrites() 





"""
Name: proc_plotPotentials_pyramidalAxon()
Type: Procedure

Input Parameter(s): 1. Object - matplotlib pyplot (plt)
                    2. List of recording locations (segments) along the axon
                    3. List of membrane potential vectors for the axon
                    4. Time vector
                    5. Start location (0 or 1)
                    6. Stop location (1 or 2)
Return Parameter: None

Description: Procedure to plot membrane potential(s) along the axon.

Edit History: Created by Nilapratim Sengupta in August-September 2023.
"""
def proc_plotPotentials_pyramidalAxon(plt, list_recLocations_axon, vecList_potentials_axon, vec_time, location_start_axon, location_stop_axon):  
    
    # Creating trace, label and indices lists
    list_traces_start = []
    list_traces_stop = []
    list_labels_start = []
    list_labels_stop = []
    list_indices_proximal = [1, 2, 5, 8, 11]
    list_indices_medial = [1, 3, 6, 9, 12]
    list_indices_distal = [1, 4, 7, 10, 13]
    
    # Selecting membrane potential vectors for start location
    if (location_start_axon == 0):
        for index in list_indices_proximal:
            list_traces_start.append(vecList_potentials_axon[index])
            list_labels_start.append(list_recLocations_axon[index])

    elif (location_start_axon == 1):
        for index in list_indices_medial:
            list_traces_start.append(vecList_potentials_axon[index])
            list_labels_start.append(list_recLocations_axon[index])
        
    # Selecting membrane potential vectors for stop location
    if (location_stop_axon == 1):
        for index in list_indices_medial:
            list_traces_stop.append(vecList_potentials_axon[index])
            list_labels_stop.append(list_recLocations_axon[index])
            
    elif (location_stop_axon == 2):
        for index in list_indices_distal:
            list_traces_stop.append(vecList_potentials_axon[index])
            list_labels_stop.append(list_recLocations_axon[index])

    # Creating a figure
    fig0, ax0 = plt.subplots(nrows=2, sharex=True, squeeze=False)
    ax0_1, ax0_2 = ax0.flatten()
    
    # Plotting in the figure
    for trace, label in zip(list_traces_start, list_labels_start):
        ax0_1.plot(vec_time, trace, label=label)
    for trace, label in zip(list_traces_stop, list_labels_stop):
        ax0_2.plot(vec_time, trace, label=label)

    # Adding figure details
    ax0_2.set_xlabel('Time (ms)')
    ax0_1.set_ylabel('Membrane Potential (mV)')
    ax0_2.set_ylabel('Membrane Potential (mV)')
    ax0_1.set_title('Axon - Location 1')
    ax0_2.set_title('Axon - Location 2')
    ax0_1.legend()
    ax0_2.legend()

    # Visualizing the figure
    # plt.show()
    
# End of proc_plotPotentials_pyramidalAxon()  





"""
Name: proc_plotPotentials_pyramidalCollateral()
Type: Procedure

Input Parameter(s): 1. Object - matplotlib pyplot (plt)
                    2. List of recording locations (segments) along the axon
                    3. List of membrane potential vectors for the axon
                    4. List of recording locations (segments) along the collateral(s)
                    5. List of membrane potential vectors for the collateral(s)
                    6. Number of collaterals
                    7. Time vector
                    8. Stop time
                    9. Start location (-1 or 0 or 1)
                   10. Stop location (1 or 2)
Return Parameter: None

Description: Procedure to plot membrane potential(s) along the collateral(s).

Edit History: Created by Nilapratim Sengupta in September 2023.
              Edited by Nilapratim Sengupta in December 2023 to update options for location_start_collateral.
"""
def proc_plotPotentials_pyramidalCollateral(plt, list_recLocations_axon, vecList_potentials_axon, list_recLocations_collaterals, vecList_potentials_collaterals, count_collaterals, vec_time, time_stop, location_start_collateral, location_stop_collateral):  
    
    # Looping over the collateral(s)
    for loopCounter_collateral in range(count_collaterals):
        
        index_collateral = loopCounter_collateral
            
        # Creating trace, label and indices lists
        list_indices_axon = [1, 2, 5, 8, 11]
        list_traces_start = []
        list_traces_stop = []
        list_labels_start = []
        list_labels_stop = []
        offset = index_collateral*12
        list_indices_proximal = [offset+0, offset+3, offset+6, offset+9]
        list_indices_medial = [offset+1, offset+4, offset+7, offset+10]
        list_indices_distal = [offset+2, offset+5, offset+8, offset+11]
        
        # Selecting membrane potential vectors for start location
        if (location_start_collateral == -1):
            for index in list_indices_axon:
                list_traces_start.append(vecList_potentials_axon[index])
                list_labels_start.append(list_recLocations_axon[index])
            
        elif (location_start_collateral == 0):
            for index in list_indices_proximal:
                list_traces_start.append(vecList_potentials_collaterals[index])
                list_labels_start.append(list_recLocations_collaterals[index])

        elif (location_start_collateral == 1):
            for index in list_indices_medial:
                list_traces_start.append(vecList_potentials_collaterals[index])
                list_labels_start.append(list_recLocations_collaterals[index])
            
        # Selecting membrane potential vectors for stop location
        if (location_stop_collateral == 1):
            for index in list_indices_medial:
                list_traces_stop.append(vecList_potentials_collaterals[index])
                list_labels_stop.append(list_recLocations_collaterals[index])
                
        elif (location_stop_collateral == 2):
            for index in list_indices_distal:
                list_traces_stop.append(vecList_potentials_collaterals[index])
                list_labels_stop.append(list_recLocations_collaterals[index])

        # Creating a figure
        fig0, ax0 = plt.subplots(nrows=2, sharex=True, squeeze=False)
        ax0_1, ax0_2 = ax0.flatten()
        
        # Plotting in the figure
        for trace, label in zip(list_traces_start, list_labels_start):
            ax0_1.plot(vec_time, trace, label=label)
        for trace, label in zip(list_traces_stop, list_labels_stop):
            ax0_2.plot(vec_time, trace, label=label)

        # Adding figure details
        ax0_2.set_xlabel('Time (ms)')
        ax0_1.set_ylabel('Membrane Potential (mV)')
        ax0_2.set_ylabel('Membrane Potential (mV)')
        ax0_1.set_title(f'Collateral {index_collateral} - Location 1')
        ax0_2.set_title(f'Collateral {index_collateral} - Location 2')
        ax0_1.legend()
        ax0_2.legend()

    # Visualizing the figure
    # plt.show()
    
# End of proc_plotPotentials_pyramidalCollateral() 





"""
Name: proc_savePotentials_pyramidalSomaDendrites()
Type: Procedure

Input Parameter(s): 1. Object - pandas dataframe (pd)
                    2. List of recording locations (segments) in the soma and dendrites
                    3. List of membrane potential vectors for the soma and dendrites
                    4. Time vector
                    5. Filename
Return Parameter: None

Description: Procedure to save membrane potential(s) in the soma and dendrite(s) into a file.

Edit History: Created by Nilapratim Sengupta in September 2023.
"""
def proc_savePotentials_pyramidalSomaDendrites(pd, list_recLocations_somaDendrites, vecList_potentials_somaDendrites, vec_time, outputFilename_somaDendritesTraces):  

    # Creating a dictionary to hold the data
    dictionary_data = {'Time': list(vec_time)}

    # Adding membrane potential data to the dictionary
    for loc, vec in zip(list_recLocations_somaDendrites, vecList_potentials_somaDendrites):
        dictionary_data[f'{loc}'] = list(vec)

    # Creating a DataFrame from the dictionary
    dataFrame = pd.DataFrame(dictionary_data)

    # Writing the DataFrame to an Excel file
    dataFrame.to_excel(outputFilename_somaDendritesTraces, index=False)
    
# End of proc_savePotentials_pyramidalSomaDendrites()





"""
Name: proc_savePotentials_pyramidalAxon()
Type: Procedure

Input Parameter(s): 1. Object - pandas dataframe (pd)
                    2. List of recording locations (segments) along the axon
                    3. List of membrane potential vectors for the axon
                    4. Time vector
                    5. Filename
Return Parameter: None

Description: Procedure to save membrane potential(s) along the axon into a file.

Edit History: Created by Nilapratim Sengupta in September 2023.
"""
def proc_savePotentials_pyramidalAxon(pd, list_recLocations_axon, vecList_potentials_axon, vec_time, outputFilename_axonTraces):  

    # Creating a dictionary to hold the data
    dictionary_data = {'Time': list(vec_time)}

    # Adding membrane potential data to the dictionary
    for loc, vec in zip(list_recLocations_axon, vecList_potentials_axon):
        dictionary_data[f'{loc}_index'] = list(vec)

    # Creating a DataFrame from the dictionary
    dataFrame = pd.DataFrame(dictionary_data)

    # Writing the DataFrame to an Excel file
    dataFrame.to_excel(outputFilename_axonTraces, index=False)
    
# End of proc_savePotentials_pyramidalAxon()





"""
Name: proc_savePotentials_pyramidalCollateral()
Type: Procedure

Input Parameter(s): 1. Object - pandas dataframe (pd)
                    2. List of recording locations (segments) along the collateral(s)
                    3. List of membrane potential vectors for the collateral(s)
                    4. Time vector
                    5. Filename
Return Parameter: None

Description: Procedure to save membrane potential(s) along the collateral(s) into a file.

Edit History: Created by Nilapratim Sengupta in September 2023.
"""
def proc_savePotentials_pyramidalCollateral(pd, list_recLocations_collaterals, vecList_potentials_collaterals, vec_time, outputFilename_collateralsTraces):  

    # Creating a dictionary to hold the data
    dictionary_data = {'Time': list(vec_time)}

    # Adding membrane potential data to the dictionary
    counter = 0
    index = 0
    for loc, vec in zip(list_recLocations_collaterals, vecList_potentials_collaterals):
        # Ensuring data does not get overwritten for two similar locations on different collaterals
        if counter >= 12: # 12 recording locations per collateral
            counter = 0
            index += 1
        dictionary_data[f'{loc}_{index}'] = list(vec) # If two locations on the same collateral coincide (due to short collateral size), the latter overwrites the former
        counter += 1

    # Creating a DataFrame from the dictionary
    dataFrame = pd.DataFrame(dictionary_data)

    # Writing the DataFrame to an Excel file
    dataFrame.to_excel(outputFilename_collateralsTraces, index=False)
    
# End of proc_savePotentials_pyramidalCollateral()





"""
Name: proc_myInit()
Type: Procedure

Input Parameter(s): 1. Object - pyramidalCell
                    2. Object - voltageStimulus
                    3. Object - holdingStimulus
                    4. Object - currentStimulus
                    5. Initial (desired) membrane potential 
                    6. Initial duration
                    7. Stop time
Return Parameter: None

Description: Procedure for custom initialization.
                 
Edit History: Created by Nilapratim Sengupta in August 2023.
"""
def proc_myInit(pyramidalCell, voltageStimulus, holdingStimulus, currentStimulus, membranePotential_initial, duration_initial, time_stop):
    
    h.v_init = membranePotential_initial
    currentStimulus.delay = 1e9
    holdingStimulus.amp = 0

    # Invoking finitialize()
    h.finitialize(membranePotential_initial)

    # Saving dt and tstop for actual simulation duration
    dt_save = h.dt
    tstop_save = h.tstop

    # Setting 't' for holding phase and enabling variable time-step computation for this phase
    h.t = -5000
    h.tstop = h.t + duration_initial
    cvode_status = h.cvode.active()
    if cvode_status != 0:
        h.cvode.active(0)

    # Configuring the point process
    voltageStimulus.rs = 0.01
    voltageStimulus.toff = 0  # Should prevent it from delivering nonzero current when t>0.
    voltageStimulus.amp = membranePotential_initial

    # First while-loop
    print('\n\nInitialization Stage I in progress...\n') # Display for user
    while h.t < h.tstop:
        h.fadvance()

    # Re-configuring the point process
    holdingStimulus.amp = voltageStimulus.i
    voltageStimulus.rs = 1e9  # To ensure the current it delivers during a run is miniscule.
    # This is a "suspenders & belt" approach because voltageStimulus.toff = 0

    # Restoring dt and enabling holding phase for some more time using fixed time-step protocol
    h.dt = dt_save
    h.tstop = 0
    h.t = -200

    # Second while-loop
    print('\nInitialization Stage II in progress...This may take a while!\n') # Display for user
    while h.t < h.tstop:
        h.fadvance()

    # Restoring t, tstop for actual simulation
    h.t = 0
    h.tstop = tstop_save

    # In accordance with Jennifer Luebke's electrophysiology protocol
    currentStimulus.delay = 15

    # Restore and re-init cvode if necessary
    if cvode_status != 0:
        h.cvode.active(1)
        h.cvode.re_init()
    else:
        h.fcurrent()
    h.frecord_init()

    print('\nSimulation begins...\n') # Display for user
    
# End of proc_myInit()   





"""
Name: proc_parallelizeRun()
Type: Procedure

Input Parameter(s): None
Return Parameter: None

Description: Procedure to parallelize simulation.
                 
Edit History: Created by Nilapratim Sengupta in October 2023.
"""
def proc_parallelizeRun():
    
    # Measuring real time taken to simulate
    realTime_start = datetime.now()
    
    # Parameters for parallel run handling
    count_jobs = 0
    jobID = 1 # Job ID starts at 1 not 0 because 0 is master job 
    
    # Other key parameters
    type_cell = fileParameters.type_cell
    index_parameter = fileParameters.index_parameter
    method_spineCorrection = fileParameters.method_spineCorrection
    
    membranePotential_initial = fileParameters.membranePotential_initial
    duration_initial = fileParameters.duration_initial

    time_currentStart = fileParameters.time_currentStart
    time_currentStop = fileParameters.time_currentStop
    time_start = time_currentStart
    time_stop = time_currentStop + fileParameters.duration_extra
    time_step = fileParameters.time_step
    location_start_axon = fileParameters.location_start_axon
    location_stop_axon = fileParameters.location_stop_axon
    location_start_collateral = fileParameters.location_start_collateral
    location_stop_collateral = fileParameters.location_stop_collateral
    
    code_perturbation = fileParameters.code_perturbation
    code_paranodeScaling = fileParameters.code_paranodeScaling
    
    count_collaterals = fileParameters.count_collaterals
    if (count_collaterals != 0):
        connectionSiteList_collaterals = fileParameters.connectionSiteList_collaterals
        scaleFactorList_collateralDiameter = fileParameters.scaleFactorList_collateralDiameter
    
    inputFilename_axonParameters = fileParameters.inputFilename_axonParameters
    inputFilename_axonPerturbations = fileParameters.inputFilename_axonPerturbations
    if (count_collaterals != 0):
        if (count_collaterals == 1):
            inputFilename_collateralPerturbations_0 = fileParameters.inputFilename_collateralPerturbations_0
        if (count_collaterals == 2): 
            inputFilename_collateralPerturbations_0 = fileParameters.inputFilename_collateralPerturbations_0   
            inputFilename_collateralPerturbations_1 = fileParameters.inputFilename_collateralPerturbations_1
     
    currentLevel_initial = fileParameters.currentLevel_initial
    currentLevel_final = fileParameters.currentLevel_final
    currentLevel_step = fileParameters.currentLevel_step

    year = datetime.now().year
    month = datetime.now().month
    day = datetime.now().day
    simulationID = fileParameters.simulationID
    
    specialDescription = fileParameters.specialDescription

    # Creating file handlers
    outputFilenameIdentifier_results = f'{year}_{month:02d}_{day:02d}_{simulationID}'
    outputFilename_axonOutput = f'{outputFilenameIdentifier_results}_axonOutput.xlsx'
    if (count_collaterals != 0):
        outputFilename_collateralsOutput = f'{outputFilenameIdentifier_results}_collateralsOutput.xlsx'
    outputFilename_details = f'{outputFilenameIdentifier_results}_detailedOutput.txt'
    
    # Printing simulation details to file
    with open(outputFilename_details, 'w') as writer_details:
        
        # Printing the header
        writer_details.write('Printing few simulation details for record keeping.\n')
        writer_details.write('===================================================\n')
        
        # Printing special description
        writer_details.write(f'{specialDescription}')    
        
        # Printing input file names
        writer_details.write(f'\n\nPopulation File = {inputFilename_axonParameters}\n')
        writer_details.write(f'Axon Perturbation File = {inputFilename_axonPerturbations}\n')
        if (count_collaterals != 0):
            if (count_collaterals == 1):
                writer_details.write(f'Collateral Perturbation File 1 = {inputFilename_collateralPerturbations_0}\n')
            if (count_collaterals == 2):
                writer_details.write(f'Collateral Perturbation File 1 = {inputFilename_collateralPerturbations_0}\n')
                writer_details.write(f'Collateral Perturbation File 2 = {inputFilename_collateralPerturbations_1}\n')


    # Obtaining axon parameters from file
    matrix_axonParameters = fileProceduresFunctions.func_readParametersFromFile(inputFilename_axonParameters)
    size_population = matrix_axonParameters.shape[0]
    count_axonParameters = matrix_axonParameters.shape[1]
    
    # Displaying on the console
    print(f'\nPopulation Size = {size_population}') # Console-display
    print(f'No. of Parameters = {count_axonParameters}') # Console-display

    # Setting up each MPI run 
    parallelContext.runworker()
    
    # Iterating over all member(s) of the population to obtain individual axon parameters
    for loopCounter_population in range(len(matrix_axonParameters)):
        
        id_axon = matrix_axonParameters[loopCounter_population][0]
        diameter_axon = matrix_axonParameters[loopCounter_population][1]
        length_node = matrix_axonParameters[loopCounter_population][2]
        length_segment = matrix_axonParameters[loopCounter_population][3]
        count_lamellae = matrix_axonParameters[loopCounter_population][4]
        thickness_lamella = matrix_axonParameters[loopCounter_population][5]
        scaleFactor_gPas = matrix_axonParameters[loopCounter_population][6]
        scaleFactor_gNa = matrix_axonParameters[loopCounter_population][7]
        scaleFactor_gK = matrix_axonParameters[loopCounter_population][8]
        
        # Displaying on the console
        print(f'\nPopulation: {loopCounter_population}') # Console-display
        print('==============\n') # Console-display
        print(f'Axon ID = {id_axon}') # Console-display
        print(f'Axon Diameter = {diameter_axon}') # Console-display
        print(f'Node Length = {length_node}') # Console-display
        print(f'Myelinated Segment Length = {length_segment}') # Console-display
        print(f'No. of Lamellae = {count_lamellae}') # Console-display
        print(f'Thickness of each Lamella = {thickness_lamella}') # Console-display
        print(f'Scale Factor for gPas = {scaleFactor_gPas}') # Console-display
        print(f'Scale Factor for gNa = {scaleFactor_gNa}') # Console-display
        print(f'Scale Factor for gK = {scaleFactor_gK}') # Console-display
        
        
        # Obtaining axon perturbations from file
        matrix_axonPerturbations = fileProceduresFunctions.func_readPerturbationsFromFile(inputFilename_axonPerturbations)
        count_axonPerturbationColumns = matrix_axonPerturbations.shape[1]
        count_axonPerturbations = (count_axonPerturbationColumns // 2) # Integer division

        # Displaying on the console
        print(f'\nNo. of Perturbation Schemes for the Axon = {count_axonPerturbations}') # Console-display
            

        # Iterating over all perturbation scheme(s)
        for loopCounter_perturbation in range (count_axonPerturbations): # Ensure that the axon and the collateral(s) have the same number of perturbations
            

            # Iterating over all current level(s)
            loopCounter_currentLevel = currentLevel_initial
            while loopCounter_currentLevel <= currentLevel_final:
                
                # Display for user
                print(f'\nJob ID: {jobID}\tPrepping for Population: {loopCounter_population}\tPerturbation: {loopCounter_perturbation}\tInjection: {loopCounter_currentLevel} nA.') 
                
                # Creating a list of input arguments
                list_arguments = [jobID, loopCounter_population, loopCounter_perturbation, loopCounter_currentLevel]
                
                # Invoking the function to run the simulation
                parallelContext.submit(fileProceduresFunctions.func_runSimulation, list_arguments) # Parallel mode
                # list_returnObjects = fileProceduresFunctions.func_runSimulation(list_arguments) # Serial mode
                
                # Updating job counter
                count_jobs += 1
                jobID += 1
                
                # Incrementing the currentLevel loopCounter for the while loop   
                loopCounter_currentLevel += currentLevel_step
    
    # Creating lists to hold results
    list_result_somaDendrites = [[] for _ in range(count_jobs)]
    list_result_axon = [[] for _ in range(count_jobs)]
    list_result_collaterals = [[] for _ in range(count_collaterals)]
    for loopCounter_collateral in range(count_collaterals):
        list_result_collaterals[loopCounter_collateral] = [[] for _ in range(count_jobs)]
    

    # Obtaining results from MPI jobs
    while parallelContext.working():
        
        # Obtaining the input arguments for the job
        list_arguments = parallelContext.upkpyobj()
        jobID = list_arguments[0]
        loopCounter_population = list_arguments[1]
        loopCounter_perturbation = list_arguments[2]
        loopCounter_currentLevel = list_arguments[3]
        
        # Display for the user
        print(f'\nJob ID: {jobID}\tRetrieved for Population: {loopCounter_population}\tPerturbation: {loopCounter_perturbation}\tInjection: {loopCounter_currentLevel} nA.') # Display for the user
        
        # Receiving data based on specific jobID
        parallelContext.take(jobID) 
        list_returnObjects = parallelContext.upkpyobj()
        
        # Appending the received data pertaining to the job ID to the result lists
        list_result_somaDendrites[jobID-1].append(list_returnObjects[0]) # populationID
        list_result_somaDendrites[jobID-1].append(list_returnObjects[1]) # axonID
        list_result_somaDendrites[jobID-1].append(list_returnObjects[2]) # perturbationID
        list_result_somaDendrites[jobID-1].append(list_returnObjects[3]) # currentLevel
        list_result_somaDendrites[jobID-1].append(list_returnObjects[4][0])
        list_result_somaDendrites[jobID-1].append(list_returnObjects[4][1])
        list_result_somaDendrites[jobID-1].append(list_returnObjects[4][2])
        list_result_somaDendrites[jobID-1].append(list_returnObjects[4][3])
        list_result_somaDendrites[jobID-1].append(list_returnObjects[4][4])
        list_result_somaDendrites[jobID-1].append(list_returnObjects[4][5])
        list_result_somaDendrites[jobID-1].append(list_returnObjects[4][6])

        list_result_axon[jobID-1].append(list_returnObjects[0]) # populationID
        list_result_axon[jobID-1].append(list_returnObjects[1]) # axonID
        list_result_axon[jobID-1].append(list_returnObjects[2]) # perturbationID
        list_result_axon[jobID-1].append(list_returnObjects[3]) # currentLevel
        list_result_axon[jobID-1].append(list_returnObjects[5][0])
        list_result_axon[jobID-1].append(list_returnObjects[5][1])
        list_result_axon[jobID-1].append(list_returnObjects[5][2])
        list_result_axon[jobID-1].append(list_returnObjects[5][3])
        list_result_axon[jobID-1].append(list_returnObjects[5][4])
        list_result_axon[jobID-1].append(list_returnObjects[5][5])
        list_result_axon[jobID-1].append(list_returnObjects[5][6])
        
        if (count_collaterals != 0):
            for loopCounter_collateral in range(count_collaterals):
                list_result_collaterals[loopCounter_collateral][jobID-1].append(list_returnObjects[0]) # populationID
                list_result_collaterals[loopCounter_collateral][jobID-1].append(list_returnObjects[1]) # axonID
                list_result_collaterals[loopCounter_collateral][jobID-1].append(loopCounter_collateral) # collateralID
                list_result_collaterals[loopCounter_collateral][jobID-1].append(list_returnObjects[2]) # perturbationID
                list_result_collaterals[loopCounter_collateral][jobID-1].append(list_returnObjects[3]) # currentLevel
                list_result_collaterals[loopCounter_collateral][jobID-1].append(list_returnObjects[6][loopCounter_collateral][0])
                list_result_collaterals[loopCounter_collateral][jobID-1].append(list_returnObjects[6][loopCounter_collateral][1])
                list_result_collaterals[loopCounter_collateral][jobID-1].append(list_returnObjects[6][loopCounter_collateral][2])
                list_result_collaterals[loopCounter_collateral][jobID-1].append(list_returnObjects[6][loopCounter_collateral][3])
                list_result_collaterals[loopCounter_collateral][jobID-1].append(list_returnObjects[6][loopCounter_collateral][4])
                list_result_collaterals[loopCounter_collateral][jobID-1].append(list_returnObjects[6][loopCounter_collateral][5])
                list_result_collaterals[loopCounter_collateral][jobID-1].append(list_returnObjects[6][loopCounter_collateral][6])
           
    # The while loop ends
    parallelContext.done()  
    
    # Printing simulation details to file
    with open(outputFilename_details, 'a') as writer_details:
            
        # Printing parameters regarding population size
        writer_details.write(f'\n\nPopulation Size = {size_population}\n')
        writer_details.write(f'No. of Axon Parameters = {count_axonParameters}\n')

        # Printing parameters regarding perturbations
        writer_details.write(f'No. of Perturbation Variations = {count_axonPerturbations}\n')
        
        # Printing collateral counts, connection sites and scale factor for diameter
        writer_details.write(f'\n\nNo. of Collaterals = {count_collaterals}\n')
        if (count_collaterals != 0):
            writer_details.write(f'Default Connection Sites for Collaterals = {connectionSiteList_collaterals} (Note: For remyelination, this may change!)\n')
            writer_details.write(f'Scale Factors for Collateral Diameters = {scaleFactorList_collateralDiameter}\n')

        # Printing cell type and other model details
        writer_details.write(f'\n\nCell Type = {type_cell} (2: Young; 3: Aged)\n')
        writer_details.write(f'Parameter Index = {index_parameter}\n')
        writer_details.write(f'Spine Correction Method = {method_spineCorrection}\n')
        writer_details.write(f'Code for Paranode Scaling = {code_paranodeScaling}\n')
        
        # Printing other simulation details
        writer_details.write(f'\n\nInitial Membrane Potential = {membranePotential_initial}\n')
        writer_details.write(f'Initial Duration = {duration_initial}\n')
        writer_details.write(f'\n\nTime Start = {time_start}\n') 
        writer_details.write(f'Time Stop = {time_stop}\n') # 100 ms more than time_currentStop
        writer_details.write(f'Time Step = {time_step}\n')
        writer_details.write(f'\n\nCurrent Level Start = {currentLevel_initial}\n')
        writer_details.write(f'Current Level Stop = {currentLevel_final}\n')
        writer_details.write(f'Current Level Step = {currentLevel_step}\n')
        writer_details.write(f'\n\nLocation Start Axon = {location_start_axon} (Either 0 or 1)\n')
        writer_details.write(f'Location Stop Axon = {location_stop_axon} (Either 1 or 2)\n')
        writer_details.write(f'\n\nLocation Start Collateral = {location_start_collateral} (Either -1 or 0 or 1)\n')
        writer_details.write(f'Location Stop Collateral = {location_stop_collateral} (Either 1 or 2)\n')

   
    # Specifying column headings for the axon file
    columnHeadings_axon = ['populationID', 'axonID', 'perturbationID', 'currentLevel', 'startFR', 'stopFR', 'startAPC', 'stopAPC', 'APF', 'CV', 'outcomeFlag']
    
    # Creating DataFrame from the list and writing to an Excel file
    dataFrame_axon = pd.DataFrame(list_result_axon)
    dataFrame_axon.columns = columnHeadings_axon
    with pd.ExcelWriter(outputFilename_axonOutput, engine='xlsxwriter') as writer_axon:
        dataFrame_axon.to_excel(writer_axon, sheet_name='axon', index=False)
        

    # For the collateral(s), if any   
    if (count_collaterals != 0):   
    
        # Specifying column headings for the collateral file
        columnHeadings_collaterals = ['populationID', 'axonID', 'collateralID', 'perturbationID', 'currentLevel', 'startFR', 'stopFR', 'startAPC', 'stopAPC', 'APF', 'CV', 'outcomeFlag']

        # Creating DataFrames from the dictionaries and writing to an Excel file
        with pd.ExcelWriter(outputFilename_collateralsOutput, engine='xlsxwriter') as writer_collateral:
            for loopCounter_collateral in range(count_collaterals):
                dataFrame_collateral = pd.DataFrame(list_result_collaterals[loopCounter_collateral])
                dataFrame_collateral.columns = columnHeadings_collaterals
                sheet_name = f'collateral_{loopCounter_collateral}'
                dataFrame_collateral.to_excel(writer_collateral, sheet_name=sheet_name, index=False)
    
    # Measuring real time taken to simulate
    realTime_stop = datetime.now()
    realTime_duration = (realTime_stop - realTime_start).total_seconds()
    realTime_hours, realTime_remainder = divmod(realTime_duration, 3600)
    realTime_minutes, realTime_seconds = divmod(realTime_remainder, 60)
    
    # Printing real time taken to simulate to file
    with open(outputFilename_details, 'a') as writer_details:
        writer_details.write(f'\n\nReal time taken = {realTime_hours} hrs, {realTime_minutes} mins & {realTime_seconds} secs.\n')
                
    # Displaying on the console
    print(f'\n\nSimulation {simulationID} completed!\n\n')
    
    # Exiting the program
    h.quit()
        
# End of proc_parallelizeRun()





"""
Name: func_runSimulation()
Type: Function

Input Parameter(s): List of arguments
Return Parameter: List of objects holding simulation results

Description: Function to run simulation and return results.
                 
Edit History: Created by Nilapratim Sengupta in October 2023.
              Edited by Nilapratim Sengupta in September 2024 to ensure collaterals can be modified individually.
"""
def func_runSimulation(list_arguments):
    
    # Unpacking arguments from the list
    jobID = list_arguments[0]
    loopCounter_population = list_arguments[1]
    loopCounter_perturbation = list_arguments[2]
    loopCounter_currentLevel = list_arguments[3]
    
    # Key parameters
    type_cell = fileParameters.type_cell
    index_parameter = fileParameters.index_parameter
    method_spineCorrection = fileParameters.method_spineCorrection
    
    membranePotential_initial = fileParameters.membranePotential_initial
    duration_initial = fileParameters.duration_initial

    time_currentStart = fileParameters.time_currentStart
    time_currentStop = fileParameters.time_currentStop
    time_start = time_currentStart
    time_stop = time_currentStop + fileParameters.duration_extra
    time_step = fileParameters.time_step
    location_start_axon = fileParameters.location_start_axon
    location_stop_axon = fileParameters.location_stop_axon
    location_start_collateral = fileParameters.location_start_collateral
    location_stop_collateral = fileParameters.location_stop_collateral
    
    code_perturbation = fileParameters.code_perturbation
    code_paranodeScaling = fileParameters.code_paranodeScaling
        
    count_collaterals = fileParameters.count_collaterals
    if (count_collaterals != 0):
        connectionSiteList_collaterals = fileParameters.connectionSiteList_collaterals
        scaleFactorList_collateralDiameter = fileParameters.scaleFactorList_collateralDiameter
    
    inputFilename_axonParameters = fileParameters.inputFilename_axonParameters
    inputFilename_axonPerturbations = fileParameters.inputFilename_axonPerturbations
    if (count_collaterals != 0):
        if (count_collaterals == 1):
            inputFilename_collateralPerturbations_0 = fileParameters.inputFilename_collateralPerturbations_0
        if (count_collaterals == 2): 
            inputFilename_collateralPerturbations_0 = fileParameters.inputFilename_collateralPerturbations_0   
            inputFilename_collateralPerturbations_1 = fileParameters.inputFilename_collateralPerturbations_1
    
    year = datetime.now().year
    month = datetime.now().month
    day = datetime.now().day
    simulationID = fileParameters.simulationID

    outputFilenameIdentifier_traces = f'{year}_{month:02d}_{day:02d}_{simulationID}_{loopCounter_population}_{loopCounter_perturbation}_{loopCounter_currentLevel}'

    
    # Obtaining axon parameters
    matrix_axonParameters = fileProceduresFunctions.func_readParametersFromFile(inputFilename_axonParameters)
    id_axon = matrix_axonParameters[loopCounter_population][0]
    diameter_axon = matrix_axonParameters[loopCounter_population][1]
    length_node = matrix_axonParameters[loopCounter_population][2]
    length_segment = matrix_axonParameters[loopCounter_population][3]
    count_lamellae = matrix_axonParameters[loopCounter_population][4]
    thickness_lamella = matrix_axonParameters[loopCounter_population][5]
    scaleFactor_gPas = matrix_axonParameters[loopCounter_population][6]
    scaleFactor_gNa = matrix_axonParameters[loopCounter_population][7]
    scaleFactor_gK = matrix_axonParameters[loopCounter_population][8]
    
    # Updating lists of perturbation codes and lamellae fractions for the axon
    matrix_axonPerturbations = fileProceduresFunctions.func_readPerturbationsFromFile(inputFilename_axonPerturbations)
    initialList_perturbationCode_axon = [row[2*loopCounter_perturbation] for row in matrix_axonPerturbations]
    initialList_lamellaeFraction_axon = [row[2*loopCounter_perturbation+1] for row in matrix_axonPerturbations]
    updatedList_perturbationCode_axon = [] # Taking into account new remyelinated segments, if any
    updatedList_lamellaeFraction_axon = [] # Taking into account new remyelinated segments, if any
    
    count_axonSegments = 0
    for perturbationCode_axon, lamellaeFraction_axon in zip(initialList_perturbationCode_axon, initialList_lamellaeFraction_axon):
        if perturbationCode_axon <= 1:
            count_axonSegments += 1
            updatedList_perturbationCode_axon.append(int(perturbationCode_axon))
            updatedList_lamellaeFraction_axon.append(lamellaeFraction_axon)
        else:
            count_axonSegments += int(perturbationCode_axon)
            updatedList_perturbationCode_axon.extend([int(perturbationCode_axon)] * int(perturbationCode_axon))
            updatedList_lamellaeFraction_axon.extend([lamellaeFraction_axon] * int(perturbationCode_axon))

    # Displaying on the console
    print(f'No. of Segments in the Axon = {count_axonSegments}') # Console-display
    print(f'Perturbation Code Sequence in Axon = {updatedList_perturbationCode_axon}') # Console-display
    print(f'Lamellae Fraction Sequence in Axon = {updatedList_lamellaeFraction_axon}') # Console-display
    
    # Obtaining collateral perturbations from file(s)
    if (count_collaterals != 0):
        countList_collateralSegments = []
        listInitialList_perturbationCode_collaterals = [] # List of lists, defined differently as it gets updated differently
        listInitialList_lamellaeFraction_collaterals = [] # List of lists, defined differently as it gets updated differently
        listUpdatedList_perturbationCode_collaterals = [[] for loopCounter_collateral in range(count_collaterals)] # List of lists
        listUpdatedList_lamellaeFraction_collaterals = [[] for loopCounter_collateral in range(count_collaterals)] # List of lists
        for loopCounter_collateral in range(count_collaterals):
            tempFilename_collateralPerturbations = (f'inputFilename_collateralPerturbations_{loopCounter_collateral}')
            matrix_collateralPerturbations = fileProceduresFunctions.func_readPerturbationsFromFile(locals().get(tempFilename_collateralPerturbations)) 
            listInitialList_perturbationCode_collaterals.append([row[2*loopCounter_perturbation] for row in matrix_collateralPerturbations])
            listInitialList_lamellaeFraction_collaterals.append([row[2*loopCounter_perturbation+1] for row in matrix_collateralPerturbations])
            
            # Updating lists of perturbation codes and lamellae fractions for the collateral(s)
            count_collateralSegments = 0
            for perturbationCode_collateral, lamellaeFraction_collateral in zip(listInitialList_perturbationCode_collaterals[loopCounter_collateral], listInitialList_lamellaeFraction_collaterals[loopCounter_collateral]):
                if perturbationCode_collateral <= 1:
                    count_collateralSegments += 1
                    listUpdatedList_perturbationCode_collaterals[loopCounter_collateral].append(int(perturbationCode_collateral))
                    listUpdatedList_lamellaeFraction_collaterals[loopCounter_collateral].append(lamellaeFraction_collateral)
                else:
                    count_collateralSegments += int(perturbationCode_collateral)
                    listUpdatedList_perturbationCode_collaterals[loopCounter_collateral].extend([int(perturbationCode_collateral)] * int(perturbationCode_collateral))
                    listUpdatedList_lamellaeFraction_collaterals[loopCounter_collateral].extend([lamellaeFraction_collateral] * int(perturbationCode_collateral))
        
            countList_collateralSegments.append(count_collateralSegments)

            # Displaying on the console
            print(f'\nNo. of Segments in the Collateral {loopCounter_collateral} = {countList_collateralSegments[loopCounter_collateral]}') # Console-display
            print(f'Perturbation Code Sequence in Collateral {loopCounter_collateral} = {listUpdatedList_perturbationCode_collaterals[loopCounter_collateral]}') # Console-display
            print(f'Lamellae Fraction Sequence in Collateral {loopCounter_collateral} = {listUpdatedList_lamellaeFraction_collaterals[loopCounter_collateral]}') # Console-display
            
        # Updating connection site(s) for collateral(s) in case of perturbation (remyelination)
        connectionSiteList_collaterals_updated = fileProceduresFunctions.func_updateConnectionSiteList_pyramidalCollateral(count_collaterals, connectionSiteList_collaterals, initialList_perturbationCode_axon)  

        # Updating collateral diameter
        diameterList_collateral = [temp * diameter_axon for temp in scaleFactorList_collateralDiameter]
        
    # Creating a pyramidal cell (with defined topology)
    if (count_collaterals == 0):
        pyramidalCell = fileCellCreation.PyramidalCellTemplateWithoutCollateral(count_axonSegments)
    else:
        pyramidalCell = fileCellCreation.PyramidalCellTemplate(count_axonSegments, count_collaterals, countList_collateralSegments, connectionSiteList_collaterals_updated)  
     
    # Setting the geometry (length, diameter, nseg) and biophysics for all the soma and dendrite sections
    fileProceduresFunctions.proc_setGeometry_pyramidalSomaDendrites(pyramidalCell, type_cell)
    fileProceduresFunctions.proc_setBiophysics_pyramidalSomaDendrites(pyramidalCell)
        
    # Applying spine-correction to the dendritic sections
    fileProceduresFunctions.proc_applySpineCorrection(pyramidalCell, method_spineCorrection)
    
    
    # Computing extracellular parameter values for unperturbed myelinated segments along the axon
    list_extracellularParameters_axon = fileProceduresFunctions.func_calculateExtracellularParameters(diameter_axon, length_node, length_segment, count_lamellae, thickness_lamella, code_perturbation, code_paranodeScaling)
   
    # Setting the geometry (length, diameter, nseg) and biophysics for all the sections of the axon
    fileProceduresFunctions.proc_setGeometry_pyramidalAxon(pyramidalCell, diameter_axon, length_node, length_segment, code_perturbation, code_paranodeScaling)
    fileProceduresFunctions.proc_setBiophysics_pyramidalAxon(pyramidalCell, list_extracellularParameters_axon)

    # Implementing perturbation(s), if any, along the axon
    for index_segment, (perturbationCode_axon, lamellaeFraction_axon) in enumerate(zip(updatedList_perturbationCode_axon, updatedList_lamellaeFraction_axon)):
        if (perturbationCode_axon != 0 or lamellaeFraction_axon != 1):
            fileProceduresFunctions.proc_applyPerturbation_pyramidalAxon(pyramidalCell.pyramidalAxon, index_segment, diameter_axon, length_node, length_segment, lamellaeFraction_axon*count_lamellae, thickness_lamella, perturbationCode_axon, code_paranodeScaling)

    
    # Computing extracellular parameter values for unperturbed myelinated segments along the collateral(s)
    # Setting the geometry (length, diameter, nseg) and biophysics for all the sections of the collateral(s)
    if (count_collaterals != 0):
        for loopCounter_collateral in range(count_collaterals):
            diameter_collateral = diameterList_collateral[loopCounter_collateral]
            # print(f'\ndiameter of Collateral {loopCounter_collateral} = {diameter_collateral}') # Console-display
            list_extracellularParameters_collateral = fileProceduresFunctions.func_calculateExtracellularParameters(diameter_collateral, length_node, length_segment, count_lamellae, thickness_lamella, code_perturbation, code_paranodeScaling)
            fileProceduresFunctions.proc_setGeometry_pyramidalCollateral(pyramidalCell, loopCounter_collateral, diameter_collateral, length_node, length_segment, code_perturbation, code_paranodeScaling) 
            fileProceduresFunctions.proc_setBiophysics_pyramidalCollateral(pyramidalCell, loopCounter_collateral, list_extracellularParameters_collateral)
    

    # Implementing perturbation(s), if any, along the collateral(s)
    if (count_collaterals != 0):
        for loopCounter_collateral in range(count_collaterals):
            diameter_collateral = diameterList_collateral[loopCounter_collateral]
            # print(f'\ndiameter of Collateral {loopCounter_collateral} = {diameter_collateral}') # Console-display
            for index_segment, (perturbationCode_collateral, lamellaeFraction_collateral) in enumerate(zip(listUpdatedList_perturbationCode_collaterals[loopCounter_collateral], listUpdatedList_lamellaeFraction_collaterals[loopCounter_collateral])):
                if (perturbationCode_collateral != 0 or lamellaeFraction_collateral != 1):
                    fileProceduresFunctions.proc_applyPerturbation_pyramidalCollateral(pyramidalCell.pyramidalCollateral[loopCounter_collateral], index_segment, diameter_collateral, length_node, length_segment, lamellaeFraction_collateral*count_lamellae, thickness_lamella, perturbationCode_collateral, code_paranodeScaling)

    # Setting channel conductances and kinetics
    fileProceduresFunctions.proc_setChannelConductances_pyramidalSomaDendrites(pyramidalCell, type_cell, index_parameter)
    fileProceduresFunctions.proc_setChannelConductances_pyramidalAxon(pyramidalCell, type_cell, index_parameter)
    if (count_collaterals != 0):
        fileProceduresFunctions.proc_setChannelConductances_pyramidalCollateral(pyramidalCell, count_collaterals, type_cell, index_parameter)
    fileProceduresFunctions.proc_setChannelKinetics(pyramidalCell, type_cell, index_parameter)
    fileProceduresFunctions.proc_setParameters_final(pyramidalCell, count_collaterals, scaleFactor_gPas, scaleFactor_gNa, scaleFactor_gK)
    
    # Applying current clamp
    [holdingStimulus, currentStimulus] = fileProceduresFunctions.func_setupIClamp(pyramidalCell, time_currentStart, time_currentStop)
    currentStimulus.amp = loopCounter_currentLevel
    currentStimulus.dur = time_currentStop - time_currentStart  
    currentStimulus.delay = time_currentStart

    
    # Applying voltage clamp
    voltageStimulus = fileProceduresFunctions.func_setupVClamp(pyramidalCell, membranePotential_initial)
    
    # Setting up recording locations
    list_recLocations_somaDendrites = fileProceduresFunctions.func_setRecLocations_pyramidalSomaDendrites(pyramidalCell)
    list_recLocations_axon = fileProceduresFunctions.func_setRecLocations_pyramidalAxon(pyramidalCell, count_axonSegments)
    if (count_collaterals != 0):
        list_recLocations_collaterals = fileProceduresFunctions.func_setRecLocations_pyramidalCollateral(pyramidalCell, count_collaterals, countList_collateralSegments)
    
    # Setting up vectors to record spike times
    vecList_spikeTimes_somaDendrites = fileProceduresFunctions.func_detectSpikeTimes_pyramidalSomaDendrites(list_recLocations_somaDendrites)
    vecList_spikeTimes_axon = fileProceduresFunctions.func_detectSpikeTimes_pyramidalAxon(list_recLocations_axon)
    if (count_collaterals != 0):
        vecList_spikeTimes_collaterals = fileProceduresFunctions.func_detectSpikeTimes_pyramidalCollateral(list_recLocations_collaterals)
    
    # Setting up vectors to record time and membrane potentials
    vec_time = h.Vector().record(h._ref_t)
    vecList_potentials_somaDendrites = fileProceduresFunctions.func_recordVectors_pyramidalSomaDendrites(list_recLocations_somaDendrites)
    vecList_potentials_axon = fileProceduresFunctions.func_recordVectors_pyramidalAxon(list_recLocations_axon)
    if (count_collaterals != 0):
        vecList_potentials_collaterals = fileProceduresFunctions.func_recordVectors_pyramidalCollateral(list_recLocations_collaterals)

    # Initializing the model
    fileProceduresFunctions.proc_myInit(pyramidalCell, voltageStimulus, holdingStimulus, currentStimulus, membranePotential_initial, duration_initial, time_stop)
    
    # Printing section details to console
    h.psection(sec=pyramidalCell.pyramidalSomaDendrites.comp[1])
    h.psection(sec=pyramidalCell.pyramidalSomaDendrites.comp[3])
    h.psection(sec=pyramidalCell.pyramidalSomaDendrites.comp[68])
    h.psection(sec=pyramidalCell.pyramidalAxon.hill_axon)
    h.psection(sec=pyramidalCell.pyramidalAxon.iseg_axon)
    h.psection(sec=pyramidalCell.pyramidalAxon.nodes_axon[0])
    h.psection(sec=pyramidalCell.pyramidalAxon.paranodeOnes_axon[0])
    h.psection(sec=pyramidalCell.pyramidalAxon.paranodeTwos_axon[0])
    h.psection(sec=pyramidalCell.pyramidalAxon.paranodeThrees_axon[0])
    h.psection(sec=pyramidalCell.pyramidalAxon.paranodeFours_axon[0])
    h.psection(sec=pyramidalCell.pyramidalAxon.juxtaparanodes_axon[0])
    h.psection(sec=pyramidalCell.pyramidalAxon.internodes_axon[0])
    for loopCounter_collateral in range(count_collaterals):
        print(f'\nFor Collateral {loopCounter_collateral}') # Console-display
        h.psection(sec=pyramidalCell.pyramidalCollateral[loopCounter_collateral].nodes_collateral[0])
        h.psection(sec=pyramidalCell.pyramidalCollateral[loopCounter_collateral].paranodeOnes_collateral[0])
        h.psection(sec=pyramidalCell.pyramidalCollateral[loopCounter_collateral].paranodeTwos_collateral[0])
        h.psection(sec=pyramidalCell.pyramidalCollateral[loopCounter_collateral].paranodeThrees_collateral[0])
        h.psection(sec=pyramidalCell.pyramidalCollateral[loopCounter_collateral].paranodeFours_collateral[0])
        h.psection(sec=pyramidalCell.pyramidalCollateral[loopCounter_collateral].juxtaparanodes_collateral[0])
        h.psection(sec=pyramidalCell.pyramidalCollateral[loopCounter_collateral].internodes_collateral[0])

    # Simulating the model
    h.continuerun(time_stop)

    # Obtaining propagation details along different paths
    list_propagationDetails_somaDendrites = fileProceduresFunctions.func_calculatePropagationDetails_pyramidalSomaDendrites(pyramidalCell, vecList_spikeTimes_somaDendrites, time_start, time_stop)
    list_propagationDetails_axon = fileProceduresFunctions.func_calculatePropagationDetails_pyramidalAxon(pyramidalCell, vecList_spikeTimes_axon, count_axonSegments, location_start_axon, location_stop_axon, time_start, time_stop)
    if (count_collaterals != 0):
        list_propagationDetails_collaterals = [] # [[] for loopCounter_collateral in range(count_collaterals)]
        for loopCounter_collateral in range(count_collaterals):
            index_collateral = loopCounter_collateral
            count_collateralSegments = countList_collateralSegments[loopCounter_collateral]
            list_propagationDetails_collaterals.append(fileProceduresFunctions.func_calculatePropagationDetails_pyramidalCollateral(pyramidalCell, vecList_spikeTimes_axon, vecList_spikeTimes_collaterals, index_collateral, count_collateralSegments, location_start_collateral, location_stop_collateral, time_start, time_stop))

    # # Plotting graphs
    # fileProceduresFunctions.proc_plotPotentials_pyramidalSomaDendrites(plt, list_recLocations_somaDendrites, vecList_potentials_somaDendrites, vec_time)
    # fileProceduresFunctions.proc_plotPotentials_pyramidalAxon(plt, list_recLocations_axon, vecList_potentials_axon, vec_time, location_start_axon, location_stop_axon)
    # if (count_collaterals != 0):
    #     fileProceduresFunctions.proc_plotPotentials_pyramidalCollateral(plt, list_recLocations_axon, vecList_potentials_axon, list_recLocations_collaterals, vecList_potentials_collaterals, count_collaterals, vec_time, time_stop, location_start_collateral, location_stop_collateral)
    # # Visualizing the graphs
    # plt.show()
    
    # # Saving trace data
    # outputFilename_somaDendritesTraces = f'{outputFilenameIdentifier_traces}_somaDendritesTraces.xlsx'
    # outputFilename_axonTraces = f'{outputFilenameIdentifier_traces}_axonTraces.xlsx'
    # if (count_collaterals != 0):
    #     outputFilename_collateralsTraces = f'{outputFilenameIdentifier_traces}_collateralsTraces.xlsx'
    # fileProceduresFunctions.proc_savePotentials_pyramidalSomaDendrites(pd, list_recLocations_somaDendrites, vecList_potentials_somaDendrites, vec_time, outputFilename_somaDendritesTraces)
    # fileProceduresFunctions.proc_savePotentials_pyramidalAxon(pd, list_recLocations_axon, vecList_potentials_axon, vec_time, outputFilename_axonTraces)
    # if (count_collaterals != 0):
    #     fileProceduresFunctions.proc_savePotentials_pyramidalCollateral(pd, list_recLocations_collaterals, vecList_potentials_collaterals, vec_time, outputFilename_collateralsTraces)
    
    # Creating list to be returned/posted by the worker node
    list_returnObjects = []
    list_returnObjects.append(loopCounter_population)
    list_returnObjects.append(id_axon)
    list_returnObjects.append(loopCounter_perturbation)
    list_returnObjects.append(loopCounter_currentLevel)
    list_returnObjects.append(list_propagationDetails_somaDendrites)
    list_returnObjects.append(list_propagationDetails_axon)
    if (count_collaterals != 0):
        list_returnObjects.append(list_propagationDetails_collaterals)
    
    # Posting the list updated with simulated data - for parallelized simulations
    parallelContext.post(jobID, list_returnObjects) 
    
    # Returning the list updated with simulated data
    return list_returnObjects
        
# End of func_runSimulation()