{ "cells": [ { "cell_type": "code", "execution_count": null, "id": "2ee1d062", "metadata": {}, "outputs": [], "source": [ "# This code is written by Nooshin Abdollahi\n", "# Information about this code:\n", "# - Motor axons are not included\n", "# - there are not transverse connections between Boundary and Boundary" ] }, { "cell_type": "code", "execution_count": 1, "id": "af4c646e", "metadata": {}, "outputs": [], "source": [ "# show the time of execution\n", "from datetime import datetime\n", "start_time = datetime.now()\n" ] }, { "cell_type": "code", "execution_count": 2, "id": "493e7e8a", "metadata": {}, "outputs": [], "source": [ "from neuron import h\n", "import netpyne \n", "from netpyne import specs, sim \n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "from typing import Tuple, List\n", "import math\n", "import sys\n", "\n", "\n", "%matplotlib inline" ] }, { "cell_type": "code", "execution_count": 3, "id": "d05a8722", "metadata": {}, "outputs": [], "source": [ "# Import nesseccery files from Matlab\n", "\n", "R = np.loadtxt(\"R.txt\") # All axons with different radius\n", "G = np.loadtxt(\"G.txt\") # Axon's groups\n", "C = np.loadtxt(\"C.txt\") # Coordinates of each axon (x,y)\n", "neighboringAxon = np.loadtxt(\"neighboringAxon.txt\")\n", "dist = np.loadtxt(\"dist.txt\") \n", "\n", "unique_radius = np.loadtxt(\"unique_radius.txt\") # including different types\n", "Number_of_nodes = np.loadtxt(\"Number_of_nodes.txt\") # Number of nodes for the specified axon total length\n", "\n", "parameters = np.loadtxt(\"parameters.txt\") \n", "\n", "# importing all the connections\n", "import scipy.io as io\n", "\n", "for i in range(len(unique_radius)):\n", " for j in range(len(unique_radius)):\n", " if j>=i:\n", " l = [i, j]\n", " z = ''.join([str(n) for n in l])\n", " Input = io.loadmat('Connect_types_{}.mat'.format(z) , squeeze_me=True) \n", " I = Input['SAVE']; \n", " locals()[\"Connect_types_\"+str(z)]=[]\n", " for v in range(len(I)):\n", " D = I[v].strip() \n", " locals()[\"Connect_types_\"+str(z)].append(D) \n", "\n", "\n", "# Boundary connections\n", "for i in range(len(unique_radius)):\n", " Input = io.loadmat('Boundary_to_{}.mat'.format(i) , squeeze_me=True) \n", " I = Input['SAVE']; \n", " locals()[\"Boundary_to_\"+str(i)]=[]\n", " for v in range(len(I)):\n", " D = I[v].strip() \n", " locals()[\"Boundary_to_\"+str(i)].append(D) \n", " \n", "\n", "\n", "#\n", "Boundary_coordinates = np.loadtxt(\"Boundary_coordinates.txt\")\n", "Boundary_neighboring = np.loadtxt(\"Boundary_neighboring.txt\")\n", "Boundary_dist = np.loadtxt(\"Boundary_dist.txt\") \n", "\n", "\n", "############## importing files related to transverse resistance (Rg) and Areas\n", "\n", "for i in range(len(unique_radius)):\n", " for j in range(len(unique_radius)):\n", " if j>=i:\n", " l = [i, j]\n", " z = ''.join([str(n) for n in l])\n", " Input = np.loadtxt('Rg_{}.txt'.format(z) ) \n", " locals()[\"Rg_\"+str(z)]=Input\n", " \n", "\n", "\n", " \n", "for i in range(len(unique_radius)):\n", " Input = np.loadtxt('Boundary_Rg_{}.txt'.format(i) ) \n", " locals()[\"Boundary_Rg_\"+str(i)]=Input\n", "\n", " \n", " \n", " \n", " \n", "for i in range(1,len(unique_radius)):\n", " for j in range(1,len(unique_radius)):\n", " if j>i:\n", " l = [i, j]\n", " z = ''.join([str(n) for n in l])\n", " Input = np.loadtxt('Areas_{}.txt'.format(z) ) \n", " locals()[\"Areas_\"+str(z)]=Input\n", " \n", " \n", " \n", " \n" ] }, { "cell_type": "code", "execution_count": 4, "id": "cf1c9f69", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\t1 \n", "\t1 \n", "\t1 \n", "\t1 \n", "\t1 \n", "\t1 \n", "\t1 \n", "\t1 \n", "\t1 \n", "\t1 \n" ] } ], "source": [ "# Network parameters\n", "netParams = specs.NetParams()\n", "\n", "netParams.sizeX=3000\n", "netParams.sizeY=3000\n", "netParams.sizeZ=3000\n", "\n", "\n", "################################# Importing Axons(including C fibers and the others) and Boundary ####################################\n", "\n", "netParams.importCellParams(\n", " cellInstance=True,\n", " label='Boundary', \n", " conds={'cellType': 'Boundary', 'cellModel': 'Boundary'},\n", " fileName='Boundarycable.hoc', \n", " cellName='Boundary', \n", " importSynMechs=True) ;\n", "\n", "\n", "## type 0 \n", "netParams.importCellParams(\n", " cellInstance=True,\n", " label='Cfiber', \n", " conds={'cellType': 'Cfiber', 'cellModel': 'Cfiber'},\n", " fileName='Cfiber.hoc', \n", " cellName='Cfiber', \n", " importSynMechs=True) ;\n", "\n", "\n", "\n", "# Myelinated axons have different types (i.e. diameters)\n", "# How many types... do I have? print(len(unique_radius)-1), -1 because the first eleman is for C fiber\n", "# each type is a specific diameter\n", "\n", "netParams.importCellParams(\n", " cellInstance=True,\n", " label='type1', \n", " conds={'cellType': 'type1', 'cellModel': 'type1'},\n", " fileName='type1.hoc', \n", " cellName='type1', \n", " importSynMechs=True) ;\n", "\n", "netParams.importCellParams(\n", " cellInstance=True,\n", " label='type2', \n", " conds={'cellType': 'type2', 'cellModel': 'type2'},\n", " fileName='type2.hoc', \n", " cellName='type2', \n", " importSynMechs=True) ;\n", "\n", "netParams.importCellParams(\n", " cellInstance=True,\n", " label='type3', \n", " conds={'cellType': 'type3', 'cellModel': 'type3'},\n", " fileName='type3.hoc', \n", " cellName='type3', \n", " importSynMechs=True) ;\n", "\n", "netParams.importCellParams(\n", " cellInstance=True,\n", " label='type4', \n", " conds={'cellType': 'type4', 'cellModel': 'type4'},\n", " fileName='type4.hoc', \n", " cellName='type4', \n", " importSynMechs=True) ;\n", "\n", "\n", "netParams.importCellParams(\n", " cellInstance=True,\n", " label='type5', \n", " conds={'cellType': 'type5', 'cellModel': 'type5'},\n", " fileName='type5.hoc', \n", " cellName='type5', \n", " importSynMechs=True) ;\n", "\n", " \n", "netParams.importCellParams(\n", " cellInstance=True,\n", " label='type6', \n", " conds={'cellType': 'type6', 'cellModel': 'type6'},\n", " fileName='type6.hoc', \n", " cellName='type6', \n", " importSynMechs=True) ; \n", " \n", " \n", "netParams.importCellParams(\n", " cellInstance=True,\n", " label='type7', \n", " conds={'cellType': 'type7', 'cellModel': 'type7'},\n", " fileName='type7.hoc', \n", " cellName='type7', \n", " importSynMechs=True) ; \n", "\n", "\n", "netParams.importCellParams(\n", " cellInstance=True,\n", " label='type8', \n", " conds={'cellType': 'type8', 'cellModel': 'type8'},\n", " fileName='type8.hoc', \n", " cellName='type8', \n", " importSynMechs=True) ;\n", "\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": 6, "id": "d5ef8f97", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "29\n" ] } ], "source": [ "###################################### Locating each axon in specific (x,y) #################################################\n", "\n", "\n", "for i in range(len(R)):\n", " x = np.where(unique_radius == R[i])\n", " \n", " if x[0]==0:\n", " netParams.popParams[\"Axon%s\" %i] = {\n", " 'cellType': 'Cfiber', \n", " 'numCells':1 , \n", " 'cellModel': 'Cfiber', \n", " 'xRange':[C[i][0], C[i][0]], \n", " 'yRange':[0, 0], \n", " 'zRange':[C[i][1], C[i][1]]} \n", "\n", "# x=1,...,x=len(unique_radius)-1 \n", " \n", " if x[0]==1:\n", " netParams.popParams[\"Axon%s\" %i] = {\n", " 'cellType': 'type1', \n", " 'numCells':1 , \n", " 'cellModel': 'type1', \n", " 'xRange':[C[i][0], C[i][0]], \n", " 'yRange':[0, 0], \n", " 'zRange':[C[i][1], C[i][1]]} \n", "\n", " if x[0]==2:\n", " netParams.popParams[\"Axon%s\" %i] = {\n", " 'cellType': 'type2', \n", " 'numCells':1 , \n", " 'cellModel': 'type2', \n", " 'xRange':[C[i][0], C[i][0]], \n", " 'yRange':[0, 0], \n", " 'zRange':[C[i][1], C[i][1]]} \n", " \n", " \n", " if x[0]==3:\n", " netParams.popParams[\"Axon%s\" %i] = {\n", " 'cellType': 'type3', \n", " 'numCells':1 , \n", " 'cellModel': 'type3', \n", " 'xRange':[C[i][0], C[i][0]], \n", " 'yRange':[0, 0], \n", " 'zRange':[C[i][1], C[i][1]]} \n", " \n", " if x[0]==4:\n", " netParams.popParams[\"Axon%s\" %i] = {\n", " 'cellType': 'type4', \n", " 'numCells':1 , \n", " 'cellModel': 'type4', \n", " 'xRange':[C[i][0], C[i][0]], \n", " 'yRange':[0, 0], \n", " 'zRange':[C[i][1], C[i][1]]} \n", " \n", " if x[0]==5:\n", " netParams.popParams[\"Axon%s\" %i] = {\n", " 'cellType': 'type5', \n", " 'numCells':1 , \n", " 'cellModel': 'type5', \n", " 'xRange':[C[i][0], C[i][0]], \n", " 'yRange':[0, 0], \n", " 'zRange':[C[i][1], C[i][1]]} \n", " \n", " if x[0]==6:\n", " netParams.popParams[\"Axon%s\" %i] = {\n", " 'cellType': 'type6', \n", " 'numCells':1 , \n", " 'cellModel': 'type6', \n", " 'xRange':[C[i][0], C[i][0]], \n", " 'yRange':[0, 0], \n", " 'zRange':[C[i][1], C[i][1]]} \n", " \n", " if x[0]==7:\n", " netParams.popParams[\"Axon%s\" %i] = {\n", " 'cellType': 'type7', \n", " 'numCells':1 , \n", " 'cellModel': 'type7', \n", " 'xRange':[C[i][0], C[i][0]], \n", " 'yRange':[0, 0], \n", " 'zRange':[C[i][1], C[i][1]]} \n", " \n", " if x[0]==8:\n", " netParams.popParams[\"Axon%s\" %i] = {\n", " 'cellType': 'type8', \n", " 'numCells':1 , \n", " 'cellModel': 'type8', \n", " 'xRange':[C[i][0], C[i][0]], \n", " 'yRange':[0, 0], \n", " 'zRange':[C[i][1], C[i][1]]} \n", " \n", "\n", "\n", " \n", " \n", " \n", " \n", "########################################### Locating Boundary Cables ########################################################\n", "\n", "\n", "for i in range(len(Boundary_coordinates)):\n", " \n", " netParams.popParams[\"Boundary%s\" %i] = {\n", " 'cellType': 'Boundary', \n", " 'numCells':1 , \n", " 'cellModel': 'Boundary', \n", " 'xRange':[Boundary_coordinates[i][0], Boundary_coordinates[i][0]], \n", " 'yRange':[0, 0], \n", " 'zRange':[Boundary_coordinates[i][1], Boundary_coordinates[i][1]]} \n", "\n", "\n", "\n", "# in Total, how many Cells does Netpyne generate? Length(R)+len(Boundary_coordinates)\n", "print(len(R)+len(Boundary_coordinates))\n", "\n", "from IPython.display import clear_output\n", "\n", "clear_output(wait=True)\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "4adc83be", "metadata": {}, "outputs": [], "source": [ "################################################### Stimulation ############################################################\n", "# Which group of axons do you want to stimulate?\n", "# Group1: motor axons Group2: C fibers Group3: Adelta Group4: Abeta\n", "\n", "\n", "netParams.stimSourceParams['Input1'] = {'type': 'IClamp', 'del': 1, 'dur': 0.1, 'amp': 0.6}\n", "\n", "\n", "\n", "# for i in range(len(R)): \n", "# if G[i]==4: # Group 4\n", "# netParams.stimTargetParams['Input1->\"Stim_%s\"' %i] = {'source': 'Input1', 'sec':'node_0', 'loc': 0.5, 'conds': {'pop':\"Axon%s\" %i}} \n", " \n", "netParams.stimTargetParams['Input1->Stim_1'] = {'source': 'Input1', 'sec':'node_0', 'loc': 0.5, 'conds': {'pop':\"Axon2\"}} \n", " \n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "90a2f08b", "metadata": {}, "outputs": [], "source": [ "simConfig = specs.SimConfig()\n", "simConfig.hParams = {'celsius': 37 }\n", "\n", "simConfig.dt = 0.025 # Internal integration timestep to use default is 0.025\n", "simConfig.duration = 6\n", "simConfig.recordStim = True\n", "simConfig.recordStep = 0.025 # Step size in ms to save data (e.g. V traces, LFP, etc) default is 0.1\n", "#simConfig.cache_efficient = True\n", "#simConfig.cvode_active = True\n", "# simConfig.cvode_atol=0.0001\n", "# simConfig.cvode_rtol=0.0001\n", "\n", "\n", "simConfig.recordTraces = {'V_node_0' :{'sec':'node_0','loc':0.5,'var':'v'}}\n", "simConfig.analysis['plotTraces'] = {'include': ['allCells']} # ['Axon0','Axon1']\n", "\n", "simConfig.analysis['plot2Dnet'] = True\n", "simConfig.analysis['plot2Dnet'] = {'include': ['allCells'], 'view': 'xz'}\n", "\n", "\n", "\n", "#simConfig.recordLFP = [[56.39,-4000,51.74]] # Determine the location of the LFP electrode\n", "\n", "\n", "\n", "\n", "\n", "\n", "sim.create(netParams, simConfig)\n", "\n", "\n" ] }, { "cell_type": "markdown", "id": "9045099d", "metadata": {}, "source": [ "### xraxial and transverese conductances" ] }, { "cell_type": "code", "execution_count": null, "id": "41af5705", "metadata": {}, "outputs": [], "source": [ "# Since by default Netpyne does not insert the parameters of the extracellular mechanism, I insert them in this section\n", "# this section includes \"longitudinal\" resistivities (i.e. xraxial)\n", "\n", "Total_Length=10000\n", "#Section_Length=100\n", "number_boundary = 355 #Total_Length/Section_Length \n", "number_boundary = int(number_boundary)\n", "\n", "number_Cfiber = 655 #Total_Length/Section_Length \n", "number_Cfiber = int(number_boundary)\n", "\n", "\n", "\n", "rhoa=0.7e6 \n", "mycm=0.1 \n", "mygm=0.001 \n", "\n", "space_p1=0.002 \n", "space_p2=0.004\n", "space_i=0.004\n", "\n", "\n", "\n", "\n", "############################# For Boundary Cables #################################################\n", "\n", "# soma section is just for LFP recording, LFP in Netpyne does not work if at least one section is not called soma \n", "\n", "\n", "for j in range(len(R),len(R)+len(Boundary_coordinates)):\n", " \n", " S = sim.net.cells[j].secs[\"soma\"][\"hObj\"] \n", " for seg in S:\n", " seg.xraxial[0] = 1e9\n", " seg.xraxial[1] = 1e9\n", " seg.xg[0] = 1e9\n", " seg.xg[1] = 1e9\n", " seg.xc[0] = 0\n", " seg.xc[1] = 0\n", "\n", "\n", " for i in range(number_boundary): \n", " S = sim.net.cells[j].secs[\"section_%s\" %i][\"hObj\"]\n", " for seg in S:\n", " seg.xraxial[0] = 1e9\n", " seg.xraxial[1] = 1e9\n", " seg.xg[0] = 1e9\n", " seg.xg[1] = 1e9\n", " seg.xc[0] = 0\n", " seg.xc[1] = 0\n", " \n", " \n", " \n", " \n", "\n", "############################# For C fibers #######################################################\n", "\n", "\n", "\n", "xr=0.0001 # starting value\n", "\n", "\n", "\n", "for j in range(len(R)):\n", " if G[j]==2: # Group2: C fibers\n", " \n", " S = sim.net.cells[j].secs[\"soma\"][\"hObj\"] \n", " for seg in S:\n", " seg.xraxial[0] = xr\n", " seg.xraxial[1] = xr\n", " seg.xg[0] = 1e10\n", " seg.xg[1] = 1e-9\n", " seg.xc[0] = 0\n", " seg.xc[1] = 0\n", " \n", " \n", " for i in range(number_Cfiber): \n", " S = sim.net.cells[j].secs[\"axon_%s\" %i][\"hObj\"]\n", " for seg in S:\n", " seg.xraxial[0] = xr\n", " seg.xraxial[1] = xr\n", " seg.xg[0] = 1e10\n", " seg.xg[1] = 1e-9\n", " seg.xc[0] = 0\n", " seg.xc[1] = 0\n", " \n", " \n", " \n", "\n", " \n", "############################## For myelinated sensory axons ##################################### \n", "\n", "\n", "\n", "\n", "\n", "for j in range(len(R)):\n", " if G[j]!=2: # if it is not a C fiber \n", " x = np.where(unique_radius == R[j]) \n", " x = int(x[0])\n", " nodes = Number_of_nodes[x]\n", " nodes=int(nodes)\n", " \n", " \n", " nl = parameters[x][4]\n", " nodeD = parameters[x][1]\n", " paraD1 = nodeD\n", " axonD = parameters[x][0]\n", " paraD2 = axonD\n", " \n", " Rpn0 = (rhoa*.01)/((math.pi)*((((nodeD/2)+space_p1)**2)-((nodeD/2)**2)))\n", " Rpn1 = (rhoa*.01)/((math.pi)*((((paraD1/2)+space_p1)**2)-((paraD1/2)**2)))\n", " Rpn2 = (rhoa*.01)/((math.pi)*((((paraD2/2)+space_p2)**2)-((paraD2/2)**2)))\n", " Rpx = (rhoa*.01)/((math.pi)*((((axonD/2)+space_i)**2)-((axonD/2)**2)))\n", " \n", " \n", " \n", "\n", " S = sim.net.cells[j].secs[\"soma\"][\"hObj\"]\n", " for seg in S:\n", " seg.xraxial[0] = Rpn1\n", " seg.xraxial[1] = xr\n", " seg.xg[0] = mygm/(nl*2)\n", " seg.xg[1] = 1e-9 # disconnect from ground\n", " seg.xc[0] = mycm/(nl*2)\n", " seg.xc[1] = 0\n", "\n", " \n", " for i in range(nodes):\n", " S = sim.net.cells[j].secs[\"node_%s\" %i][\"hObj\"]\n", " for seg in S:\n", " seg.xraxial[0] = Rpn0\n", " seg.xraxial[1] = xr\n", " seg.xg[0] = 1e10\n", " seg.xg[1] = 1e-9\n", " seg.xc[0] = 0\n", " seg.xc[1] = 0\n", "\n", "\n", " for i in range(2*nodes):\n", " S = sim.net.cells[j].secs[\"MYSA_%s\" %i][\"hObj\"]\n", " for seg in S:\n", " seg.xraxial[0] = Rpn1\n", " seg.xraxial[1] = xr\n", " seg.xg[0] = mygm/(nl*2)\n", " seg.xg[1] = 1e-9\n", " seg.xc[0] = mycm/(nl*2)\n", " seg.xc[1] = 0\n", "\n", "\n", " for i in range(2*nodes):\n", " S = sim.net.cells[j].secs[\"FLUT_%s\" %i][\"hObj\"]\n", " for seg in S:\n", " seg.xraxial[0] = Rpn2\n", " seg.xraxial[1] = xr\n", " seg.xg[0] = mygm/(nl*2)\n", " seg.xg[1] = 1e-9\n", " seg.xc[0] = mycm/(nl*2)\n", " seg.xc[1] = 0 \n", "\n", "\n", " for i in range(6*nodes):\n", " S = sim.net.cells[j].secs[\"STIN_%s\" %i][\"hObj\"]\n", " for seg in S:\n", " seg.xraxial[0] = Rpx\n", " seg.xraxial[1] = xr\n", " seg.xg[0] = mygm/(nl*2)\n", " seg.xg[1] = 1e-9\n", " seg.xc[0] = mycm/(nl*2)\n", " seg.xc[1] = 0\n", " \n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "##############################This section is about transverse connections between axons #####################################\n", "# *** If you do not want to include ephaptic interaction, do not run this section\n", "# To model ephaptic effect, \"LinearMechanism\" in NEURON is used.\n", "\n", "\n", "\n", "rho = 1211 * 10000 # ohm-micron\n", "\n", "count = 0\n", "\n", "for i in range(len(R)): \n", "\n", " \n", " for j in range(len(R)): \n", " \n", " if neighboringAxon[i][j]==1:\n", " \n", "\n", " a1 = np.where(unique_radius == R[i]) # find type of R[i]\n", " a1 = a1[0][0]\n", " a2 = np.where(unique_radius == R[j]) # find type of R[j]\n", " a2 = a2[0][0]\n", "\n", "\n", " NSEG = 0\n", "\n", "\n", "\n", " if a1==a2:\n", " SEC = locals()[\"Connect_types_\"+str(a1)+str(a1)]\n", " RG = locals()[\"Rg_\"+str(a1)+str(a1)]\n", " area = (math.pi)*2*unique_radius[a1]*(np.ones((len(RG),1))) # micron^2\n", " area = area * 1e-8 #cm^2\n", " b1=i\n", " b2=j\n", "\n", " if a1<a2:\n", " SEC = locals()[\"Connect_types_\"+str(a1)+str(a2)]\n", " RG = locals()[\"Rg_\"+str(a1)+str(a2)]\n", " b1=i\n", " b2=j\n", " if a1==0:\n", " area = (math.pi)*2*unique_radius[a2]*(np.ones((len(RG),1)))\n", " area = area * 1e-8 #cm^2\n", " b1=j\n", " b2=i\n", " \n", " else:\n", " area = (math.pi)*2*locals()[\"Areas_\"+str(a1)+str(a2)]\n", " area = area[ : , np.newaxis]\n", " area = area * 1e-8\n", " \n", " \n", "\n", " if a1>a2:\n", " SEC = locals()[\"Connect_types_\"+str(a2)+str(a1)]\n", " RG = locals()[\"Rg_\"+str(a2)+str(a1)]\n", " b1=j\n", " b2=i\n", " if a2==0:\n", " area = (math.pi)*2*unique_radius[a1]*(np.ones((len(RG),1)))\n", " area = area * 1e-8 #cm^2\n", " b1=i\n", " b2=j\n", " \n", " else:\n", " area = (math.pi)*2*locals()[\"Areas_\"+str(a2)+str(a1)]\n", " area = area[ : , np.newaxis]\n", " area = area * 1e-8\n", " \n", " \n", " \n", " \n", " \n", "\n", "\n", " locals()[\"sl\"+str(count)] = h.SectionList()\n", "\n", " for z1 in range(int(len(SEC)/2)): \n", "\n", " S = sim.net.cells[b1].secs[SEC[z1]][\"hObj\"]\n", " NSEG=NSEG+S.nseg\n", " locals()[\"sl\"+str(count)].append(S)\n", "\n", " for z2 in range(int(len(SEC)/2),int(len(SEC))):\n", "\n", " S = sim.net.cells[b2].secs[SEC[z2]][\"hObj\"]\n", " locals()[\"sl\"+str(count)].append(S) \n", " \n", " \n", "\n", " nsegs=int(NSEG)\n", "\n", " locals()[\"gmat\"+str(count)] =h.Matrix(2*nsegs, 2*nsegs, 2)\n", " locals()[\"cmat\"+str(count)] =h.Matrix(2*nsegs, 2*nsegs, 2)\n", " locals()[\"bvec\"+str(count)] =h.Vector(2*nsegs)\n", " locals()[\"xl\"+str(count)] =h.Vector(2*nsegs)\n", " locals()[\"layer\"+str(count)] =h.Vector(2*nsegs)\n", " locals()[\"layer\"+str(count)].fill(2) # connect layer 2\n", " locals()[\"e\"+str(count)] = h.Vector(2*nsegs)\n", "\n", " for z3 in range(2*nsegs):\n", " locals()[\"xl\"+str(count)][z3] = 0.5\n", " \n", " \n", " \n", " \n", " \n", " \n", "# d = dist[i][j]\n", "# rd = rho*d\n", "# s = ((unique_radius[a1]*2)+(unique_radius[a2]*2))/2\n", "# locals()[\"RG\"+str(count)] = np.array(RG)*s\n", "# locals()[\"Resistance\"+str(count)] = rd/locals()[\"RG\"+str(count)]\n", "# locals()[\"Conductance\"+str(count)]=[]\n", "# for z4 in range(len(locals()[\"Resistance\"+str(count)])):\n", "# locals()[\"Conductance\"+str(count)].append(1/(locals()[\"Resistance\"+str(count)][z4]*area[z4]))\n", " \n", " \n", " \n", "\n", "\n", " \n", "# for z5 in range(0,nsegs,1):\n", "\n", "# locals()[\"gmat\"+str(count)].setval(z5, z5, locals()[\"Conductance\"+str(count)][z5][0] )\n", "# locals()[\"gmat\"+str(count)].setval(z5, nsegs+z5, -locals()[\"Conductance\"+str(count)][z5][0])\n", "# locals()[\"gmat\"+str(count)].setval(nsegs+z5, z5, -locals()[\"Conductance\"+str(count)][z5][0])\n", "# locals()[\"gmat\"+str(count)].setval(nsegs+z5, nsegs+z5, locals()[\"Conductance\"+str(count)][z5][0])\n", "\n", "\n", "\n", "\n", "\n", " ge=10\n", " for z5 in range(0,nsegs,1):\n", " locals()[\"gmat\"+str(count)].setval(z5, z5, ge)\n", " locals()[\"gmat\"+str(count)].setval(z5, nsegs+z5, -ge)\n", " locals()[\"gmat\"+str(count)].setval(nsegs+z5, z5, -ge)\n", " locals()[\"gmat\"+str(count)].setval(nsegs+z5, nsegs+z5, ge)\n", "\n", "\n", "\n", "\n", " locals()[\"lm\"+str(count)] = h.LinearMechanism(locals()[\"cmat\"+str(count)], locals()[\"gmat\"+str(count)], locals()[\"e\"+str(count)], locals()[\"bvec\"+str(count)], locals()[\"sl\"+str(count)], locals()[\"xl\"+str(count)], locals()[\"layer\"+str(count)])\n", "\n", " count=count+1\n", " \n", " SEC.clear\n", " del RG\n", " del area\n", " \n", " \n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "############################### Transverse connections between Boundary cables and Axons ######################################\n", "\n", "\n", "rho = 1136e5 * 10000 # ohm-micron\n", "\n", "\n", "\n", "rows = len(Boundary_neighboring)\n", "\n", "for i in range(rows):\n", " \n", " for j in range(len(R)):\n", " \n", " if Boundary_neighboring[i][j]==1:\n", " \n", " NSEG = 0\n", "\n", " a2 = np.where(unique_radius == R[j]) # find type \n", " a2 = a2[0][0]\n", " \n", " Boundary_RG = locals()[\"Boundary_Rg_\"+str(a2)]\n", " area = (math.pi)*2*unique_radius[a2]*(np.ones((len(Boundary_RG),1)))\n", " area = area * 1e-8 #cm^2\n", " \n", "\n", " SEC = locals()[\"Boundary_to_\"+str(a2)]\n", "\n", "\n", " locals()[\"sl\"+str(count)] = h.SectionList()\n", "\n", " for z1 in range(int(len(SEC)/2)): \n", "\n", " S = sim.net.cells[j].secs[SEC[z1]][\"hObj\"]\n", " NSEG=NSEG+S.nseg\n", " locals()[\"sl\"+str(count)].append(S)\n", "\n", " for z2 in range(int(len(SEC)/2),int(len(SEC))):\n", "\n", " S = sim.net.cells[len(R)+i].secs[SEC[z2]][\"hObj\"]\n", " locals()[\"sl\"+str(count)].append(S) \n", "\n", "\n", "\n", "\n", " nsegs=int(NSEG)\n", "\n", " locals()[\"gmat\"+str(count)] =h.Matrix(2*nsegs, 2*nsegs, 2)\n", " locals()[\"cmat\"+str(count)] =h.Matrix(2*nsegs, 2*nsegs, 2)\n", " locals()[\"bvec\"+str(count)] =h.Vector(2*nsegs)\n", " locals()[\"xl\"+str(count)] =h.Vector(2*nsegs)\n", " locals()[\"layer\"+str(count)] =h.Vector(2*nsegs)\n", " locals()[\"layer\"+str(count)].fill(2) # connect layer 2\n", " locals()[\"e\"+str(count)] = h.Vector(2*nsegs)\n", "\n", " for z3 in range(2*nsegs):\n", " locals()[\"xl\"+str(count)][z3] = 0.5\n", "\n", "\n", " \n", "# d = Boundary_dist[i][j]\n", "# rd = rho*d\n", "# s = (unique_radius[a2]*2)/2\n", "# locals()[\"Boundary_RG\"+str(count)] = np.array(Boundary_RG)*s\n", "# locals()[\"Resistance\"+str(count)] = rd/locals()[\"Boundary_RG\"+str(count)]\n", "# locals()[\"Conductance\"+str(count)]=[]\n", "# for z4 in range(len(locals()[\"Resistance\"+str(count)])):\n", "# locals()[\"Conductance\"+str(count)].append(1/(locals()[\"Resistance\"+str(count)][z4]*area[z4]))\n", "\n", " \n", "# for z5 in range(0,nsegs,1):\n", "\n", "# locals()[\"gmat\"+str(count)].setval(z5, z5, locals()[\"Conductance\"+str(count)][z5][0])\n", "# locals()[\"gmat\"+str(count)].setval(z5, nsegs+z5, - locals()[\"Conductance\"+str(count)][z5][0])\n", "# locals()[\"gmat\"+str(count)].setval(nsegs+z5, z5, - locals()[\"Conductance\"+str(count)][z5][0])\n", "# locals()[\"gmat\"+str(count)].setval(nsegs+z5, nsegs+z5, locals()[\"Conductance\"+str(count)][z5][0])\n", " \n", " \n", "\n", "\n", "\n", " ge=1000 # S/cm2\n", " \n", " \n", " for z5 in range(0,nsegs,1):\n", "\n", " locals()[\"gmat\"+str(count)].setval(z5, z5, ge)\n", " locals()[\"gmat\"+str(count)].setval(z5, nsegs+z5, -ge)\n", " locals()[\"gmat\"+str(count)].setval(nsegs+z5, z5, -ge)\n", " locals()[\"gmat\"+str(count)].setval(nsegs+z5, nsegs+z5, ge)\n", "\n", "\n", "\n", "\n", " locals()[\"lm\"+str(count)] = h.LinearMechanism(locals()[\"cmat\"+str(count)], locals()[\"gmat\"+str(count)], locals()[\"e\"+str(count)], locals()[\"bvec\"+str(count)], locals()[\"sl\"+str(count)], locals()[\"xl\"+str(count)], locals()[\"layer\"+str(count)])\n", "\n", " count=count+1\n", " \n", " \n", " SEC.clear\n", " del Boundary_RG\n", " del area\n", " \n", " \n", " \n", " \n", " \n", "\n", " \n", " \n", " \n", " \n", "# from IPython.display import clear_output\n", "\n", "# clear_output(wait=True)\n", "\n", "\n", " \n", "#gmat0.printf() \n", "\n", "# for sec in sl0:\n", "# print(sec)" ] }, { "cell_type": "code", "execution_count": null, "id": "9c976afd", "metadata": {}, "outputs": [], "source": [ "\n" ] }, { "cell_type": "markdown", "id": "b2a6c256", "metadata": {}, "source": [ "#### Recordings" ] }, { "cell_type": "code", "execution_count": null, "id": "d1494f97", "metadata": {}, "outputs": [], "source": [ "## Recording vext\n", "\n", "\n", "# v1 = sim.net.cells[45].secs[\"node_0\"][\"hObj\"]\n", "# ap1 = h.Vector()\n", "# t = h.Vector()\n", "# ap1.record(v1(0.5)._ref_v)\n", "\n", "# t.record(h._ref_t)" ] }, { "cell_type": "code", "execution_count": null, "id": "ca5603a0", "metadata": { "scrolled": true }, "outputs": [], "source": [ "## Recording v and vext[0], Abeta\n", "\n", "\n", "\n", "for i1 in range(len(R)): \n", " if G[i1]==4: \n", " print(i1)\n", " F = np.where(unique_radius == R[i1]) \n", " #nodes = int (Number_of_nodes[F]-1)\n", " for i3 in range(int(Number_of_nodes[F])):\n", "\n", " locals()[\"Abeta_v\"+str(i1)+str(i3)] = sim.net.cells[i1].secs[\"node_%s\"%i3][\"hObj\"]\n", " locals()[\"Abeta_ap\"+str(i1)+str(i3)] = h.Vector()\n", " locals()[\"Abeta_ap\"+str(i1)+str(i3)].record(locals()[\"Abeta_v\"+str(i1)+str(i3)](0.5)._ref_v)\n", "# locals()[\"Abeta_v_ext\"+str(i1)] = sim.net.cells[i1].secs[\"node_%s\"%nodes][\"hObj\"]\n", "# locals()[\"Abeta_ap_ext\"+str(i1)] = h.Vector()\n", "# locals()[\"Abeta_ap_ext\"+str(i1)].record(locals()[\"Abeta_v_ext\"+str(i1)](0.5)._ref_vext[0])\n", " \n", " print(i3)\n", "# print(nodes)\n", " \n", "\n", " \n", " \n", "t = h.Vector()\n", "t.record(h._ref_t)" ] }, { "cell_type": "code", "execution_count": null, "id": "4b9344bb", "metadata": {}, "outputs": [], "source": [ "## Recording v and vext[0], Adelta\n", "\n", "\n", "\n", "# for i2 in range(len(R)): \n", "# if G[i2]==3: \n", "# F = np.where(unique_radius == R[i2]) \n", "# nodes = int (Number_of_nodes[F]-1)\n", "# locals()[\"Adelta_v\"+str(i2)] = sim.net.cells[i2].secs[\"node_%s\"%nodes][\"hObj\"]\n", "# locals()[\"Adelta_ap\"+str(i2)] = h.Vector()\n", "# locals()[\"Adelta_ap\"+str(i2)].record(locals()[\"Adelta_v\"+str(i2)](0.5)._ref_v)\n", "# locals()[\"Adelta_v_ext\"+str(i2)] = sim.net.cells[i2].secs[\"node_%s\"%nodes][\"hObj\"]\n", "# locals()[\"Adelta_ap_ext\"+str(i2)] = h.Vector()\n", "# locals()[\"Adelta_ap_ext\"+str(i2)].record(locals()[\"Adelta_v_ext\"+str(i2)](0.5)._ref_vext[0])\n", "# print(i2)\n", " \n", "\n", "\n" ] }, { "cell_type": "markdown", "id": "d83f15db", "metadata": {}, "source": [ "#### Simulate and Analyze" ] }, { "cell_type": "code", "execution_count": null, "id": "cd6d9f09", "metadata": { "scrolled": true }, "outputs": [], "source": [ "sim.simulate()\n", "sim.analyze()" ] }, { "cell_type": "code", "execution_count": null, "id": "ceb34061", "metadata": {}, "outputs": [], "source": [ "# plotting\n", "\n", "#sim.analysis.plotLFP( plots = ['timeSeries', 'locations'] , electrodes=[ 'all'], lineWidth=1000 , fontSize=14, saveFig=True)\n", "\n", "# from matplotlib import pyplot\n", "# %matplotlib inline\n", "# pyplot.plot(t, ap1 )\n", "# #pyplot.xlim((0, 10))\n", "# pyplot.show()\n" ] }, { "cell_type": "code", "execution_count": null, "id": "ddb4904a", "metadata": {}, "outputs": [], "source": [ "# show the execution time\n", "\n", "end_time = datetime.now()\n", "print('Duration: {}'.format(end_time - start_time))" ] }, { "cell_type": "markdown", "id": "8f3b15f1", "metadata": {}, "source": [ "#### saving the data" ] }, { "cell_type": "code", "execution_count": null, "id": "890baeb5", "metadata": {}, "outputs": [], "source": [ "## saving the data\n", "\n", "\n", "\n", "\n", "## writing\n", "\n", "\n", "import csv\n", "\n", "with open('locals()[\"Abeta_v\"+str(i1)+str(i3)].csv', 'w', newline='') as f:\n", " csv.writer(f).writerows(zip( t , Abeta_ap500 , Abeta_ap501 , Abeta_ap502 , Abeta_ap503, Abeta_ap504 , Abeta_ap505 , Abeta_ap506 , Abeta_ap507 , Abeta_ap508 , Abeta_ap509 , Abeta_ap5010 , Abeta_ap5011 ))\n" ] }, { "cell_type": "code", "execution_count": null, "id": "16d8bddc", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "id": "d0cdd15b", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "id": "9766ae7e", "metadata": {}, "outputs": [], "source": [ "# netParams.cellParams.keys()\n", "# netParams.cellParams['']['']" ] }, { "cell_type": "code", "execution_count": null, "id": "e19fa77c", "metadata": {}, "outputs": [], "source": [ "# pyplot.plot(t, ap1 )\n", "# #pyplot.xlim((0, 10))\n", "# pyplot.show()" ] }, { "cell_type": "code", "execution_count": null, "id": "94e4f559", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.8" } }, "nbformat": 4, "nbformat_minor": 5 }