{
"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
}