"""
This module contains the functions that set up a simulation
"""import os
import sys
import numpy as np
import json
from neuron import h
from datetime import datetime
import time
import electrodes
import anatomy
import analysis as als
import biophysics as bio
import workspace as ws
import tools as tls
defget_paths():
""" Get the paths to the different folders used by the code.
Note: In the future, I will prepare this function to get all the
paths in a recursive way. For now, the relevant sub-paths are hard-coded """# Dictionary with folders
ws.folders = {}
# Current Working Directory
ws.folders["cwd"] = os.getcwd()
# Get all the folders in cwd and sub-folders
_, subdirs = tls.explore_folder(ws.folders["cwd"])
# for item in os.listdir(ws.folders["cwd"]):# if os.path.isdir(item):for subdir in subdirs:
# Add this folder
path = os.path.join(ws.folders["cwd"], subdir)
ws.folders[subdir] = path
# Check if there's no more sub-folders here
_, subsubdirs = tls.explore_folder(path)
wholeend = len(subsubdirs) == 0# If it's not a 'whole end', dig onfor ssd in subsubdirs:
path = os.path.join(ws.folders["cwd"], subdir)
# Do nothing for now# Later on, I'll finish this function# Hard-coded part: sub-folders
ws.folders["data/load"] = os.path.join(ws.folders["data"], "load")
ws.folders["data/saved"] = os.path.join(ws.folders["data"], "saved")
ws.folders["data/results"] = os.path.join(ws.folders["data"], "results")
ws.folders["src/hoc"] = os.path.join(ws.folders["cwd"], "src/hoc")
defcreate_log_file():
""" Create a log file for the simulation """# Create the name
ws.logfile = os.path.join(ws.folders["data"], "log.txt")
# Create the file
f = open(ws.logfile, "w")
f.close()
defsanity_checks():
""" Perform some checks to make sure that everything is defined
correctly and do whatever is necessary if not """# Using the Resistor Network model implies...if ws.settings["nerve model"] == "resistor network":
# that there is ephaptic couplingifnot ws.EC["presence"]:
msg = "WARNING: Using the Resistor Network Model implies "\
+ "the presence of ephaptic coupling\nEphaptic "\
+ "coupling was activated"
ws.log(msg)
ws.EC["presence"] = True
ws.EC["use"] = "resistor network"# that the stimulation forcing can't be "pre-computed"if ws.settings["stimulation"]["method"] == "pre-computed":
msg = "ERROR: Using the Resistor Network Model implies "\
+ "that the stimulation method can't be pre-computed\n"\
+ "Please check your settings"
ws.log(msg)
ws.terminate()
defread_settings():
""" Read all the settings from the json files """# General settings# Open json file containing the information about the chosen axon model
file = os.path.join(ws.folders["settings"], "settings.json")
withopen(file, "r") as f:
ws.settings = json.load(f)
for k, v in ws.settings["ephaptic coupling"].items():
try:
ws.settings["ephaptic coupling"][k] = tls.str2bool(v)
except TypeError:
# It's probably a dictionaryfor kk, vv in v.items():
ws.settings["ephaptic coupling"][k][kk] = tls.str2bool(vv)
# Shortening for the settings for ephaptic coupling
ws.EC = ws.settings["ephaptic coupling"]
# String to bool for the rest of settingsfor k, v in ws.settings.items():
try:
v_is_str = "False"in v or"True"in v
except TypeError:
# Not a string, dictionary or list
v_is_str = False# If it was a string, then parse itif v_is_str:
ws.settings[k] = tls.str2bool(v)
# Axon model# Open json file containing the information about the chosen axon model
axon_files = {}
ws.axonmodel_settings = {}
for m in ws.settings['axon models']:
af = m + '.json'
axon_files[m] = af
withopen(os.path.join(ws.folders["models"], af), "r") as f:
ws.axonmodel_settings[m] = json.load(f)
################################################################################ PHYSICAL PROPERTIES AND GEOMETRY# Anatomy# Open json file containing the information about anatomy
file = os.path.join(ws.folders["settings"], "anatomy.json")
withopen(file, "r") as f:
ws.anatomy_settings = json.load(f)
# Electrodes
file = os.path.join(ws.folders["settings"], "electrodes.json")
withopen(file, "r") as f:
ws.electrodes_settings = json.load(f)
# Container
file = os.path.join(ws.folders["settings"], "container.json")
withopen(file, "r") as f:
ws.container_settings = json.load(f)
# Tissues
file = os.path.join(ws.folders["settings"], "tissues.json")
withopen(file, "r") as f:
ws.tissues_settings = json.load(f)
# Check that everything is fine
sanity_checks()
defphysics():
""" Set up some physical properties of the system from the settings """# Length (cm)
ws.length = ws.anatomy_settings['length']
# Internodal length to fiberD ratio# This will rarely be used in some models
ws.ilD = ws.anatomy_settings['axons']['myelinated']['internodal length to fiberD ratio']
# Temperature
h.celsius = ws.settings['temperature']
# Resistivities and thicknesses# Units:# Resistivity: Ohm*cm# Distance units: um# Resistivities
ws.rho = {}
# Thicknesses
ws.th = {}
# Nominal thickness of the perineurium (um)
ws.th['perineurium'] = ws.tissues_settings['thickness']['perineurium']
# Resistivities of the tissues. They must be in MOhm*um # (convert from Ohm*cm)
ws.rho['perineurium'] = ws.tissues_settings['resistivities']['perineurium'] * 1.e-6 * 1.e4
ws.rho['epi_trans'] = ws.tissues_settings['resistivities']['epineurium']['transverse'] * 1.e-6 * 1.e4
ws.rho['epi_longt'] = ws.tissues_settings['resistivities']['epineurium']['longitudinal'] * 1.e-6 * 1.e4
ws.rho['endo_trans'] = ws.tissues_settings['resistivities']['endoneurium']['transverse'] * 1.e-6 * 1.e4
ws.rho['endo_longt'] = ws.tissues_settings['resistivities']['endoneurium']['longitudinal'] * 1.e-6 * 1.e4# Perineurium resistance (MOhm*um2)
ws.resist_perineurium = ws.rho['perineurium'] * ws.th['perineurium']
ws.log("Perineurium resistance %f MOhm*um2:"%ws.resist_perineurium)
# INFINITE RESISTIVITY OF AN IDEALLY TOTALLY RESISTIVE MEDIUM
ws.rho['SuperResistive'] = 1.e9defmiscellanea():
""" Miscellaneous variables. It's OK to have hard-coded
variables here """# DISCRETISATION OF THE VOLUME CONDUCTOR# Cables# NAELC. All this will need to be inside a json file# Number of segments for the cables
ws.nseg = ws.settings["miscellanea"]["NAELC def nseg"]
# Diameter (um)
ws.cdiam = 1.# Extracellular layers
ws.NAELC_extlayers = 1# Other variables
ws.raxmin = 1e99
ws.deltax_min = 1e99
ws.g_min = 1e99# Counters
ws.counters = {
"cables": 0,
"axons": 0,
"NAELC": 0,
"current injections": 0,
"recordings": 0
}
ws.maxnseg = 0
ws.totaldur = 0.
ws.injections = []
ws.injections_info = []
ws.recordings = []
defsetup_simcontrol():
""" Try to set up the simulation control if there's a simulation
meant to happen """try:
scsettings = ws.settings['simulation control']
h.v_init = scsettings['v_init']
ws.dt = h.dt = scsettings['dt']
ws.nt = scsettings['nt']
ws.tstop = scsettings['nt'] * h.dt
h.tstop = ws.tstop
ws.tarray = np.arange(0., h.tstop, h.dt)
except LookupError:
# There's not supposed to be any simulationpassdefprepare_workspace(remove_previous=True):
""" Prepare all the necessary information in the workspace to be
able to do the necessary work, including running a simulation """# Data management# Prepare the folders that the code will use
get_paths()
# Remove results from previous simulation(s)if remove_previous:
als.remove_previous_results()
# Create the log file
create_log_file()
# Pass a couple of functions to ws
ws.log = log
ws.terminate = terminate
# Read all the settings from the json files
read_settings()
now = datetime.now()
ws.log("Execution started on: %s"%now.strftime("%d %h %Y at %H:%M:%S"))
# Set up some physical properties of the system from the settings
physics()
# Set up a simulation
setup_simcontrol()
# Miscellaneous variables
miscellanea()
defprepare_simulation(prep_workspace=True):
""" Do the whole process of preparing a simulation """# Prepare the necessary variables from the configuration filesif prep_workspace:
prepare_workspace()
# Shortenings
stimmethod = ws.settings["stimulation"]["method"]
nervemodel = ws.settings["nerve model"]
# Create the NERVE TOPOLOGY. Cross-section only
anatomy.create_nerve()
# BIOPHYSICS# LOAD NEURON HOC CODE AND DEFINE FUNCTIONS# Load NEURON code
h.load_file(os.path.join(ws.folders["src/hoc"], "ephap.hoc"))
h.load_file(os.path.join(ws.folders["src/hoc"], "wire.hoc"))
for m in ws.axonmodel_settings.values():
h.load_file(m['hoc file'])
# BUILD ALL THE WIRES with all their length
bio.build_cables()
# Finish setting up information for the fasciclesif nervemodel == "resistor network":
bio.finish_fascicles()
# Create electrodes# This needs to be after the creation of the cables and before # connecting the ephapsesif stimmethod == "from electrodes":
ws.log("Creating the electrodes")
electrodes.create_electrodes()
# elif stimmethod == "pre-computed":# pass# EPHAPSESif ws.EC["presence"]:
if nervemodel == "resistor network":
# Create an instance of the Resistor Network
ws.log("Creating the Resistor Network")
ws.rnet = bio.ResistorNetwork()
# Get the resistances of the connections for each pair
ws.log("\tGetting resistor values")
ws.rnet.get_resistances()
# Connect resistors with ephapses in NEURON
ws.log("\tConnecting resistors")
ws.rnet.connect_resistors()
# OTHER CONNECTIONS
ws.log("Preparing connections to ground on the boundaries and xraxial")
# Prepare connections to ground in the systemif stimmethod == "from electrodes":
electrodes.prepare_ground_paths()
# Set xraxial and make the connections to groundif nervemodel == "resistor network":
bio.electrics()
# STIMULATION
ws.log("Setting up stimulation")
# Inject currents from the cuff electrode# Add injected currents, delays and durationsif stimmethod == "from electrodes":
# if nervemodel == "resistor network":
electrodes.set_stimulation()
elif stimmethod == "pre-computed":
bio.set_direct_extracellular_stimulation()
# RECORDINGS
ws.log("Setting up recordings")
# Time
als.set_time_recording()
# From electrodes# Flag: presence of electrode recordings# Set it to False by default; it will become True if a recording # electrode is created
ws.settings["electrode recordings"] = Falseif stimmethod == "from electrodes":
if nervemodel == "resistor network":
electrodes.set_recordings()
# Directly on the cables/axons
als.record_cables()
deflog(msg):
""" Write a message into the log file """withopen(ws.logfile, "a") as f:
f.write(msg + "\n")
# Print if wantedif ws.settings["verbose"]:
print(msg)
defterminate():
""" Terminate the whole program """
now = datetime.now()
log("THE PROGRAM WAS TERMINATED ON: %s"%now.strftime("%d %h %Y at %H:%M:%S"))
sys.exit()
defrun():
""" Run a simulation """if h.tstop <= ws.totaldur:
msg = 'WARNING: The simulation tstop is shorter than the total '\
+ 'stimulation time. This will cause trouble when processing '\
+ 'the recordings.'
ws.log(msg)
definitialize():
""" Initialize the simulation for the first time step """
h.finitialize(h.v_init)
h.fcurrent()
defgo_by_steps():
""" Main driver for the simulation """
initialize()
# Loop over time stepsfor i_t, t inenumerate(ws.tarray):
h.fadvance()
# Rest for 500 ms to allow the computer to cool downifFalse:
time.sleep(0.5)
# Time control
t1 = datetime.now()
msg = "Time step %i of %i. Total elapsed time: %i seconds"%(i_t, ws.nt, (t1 - t0).total_seconds())
# Log the time step into the simulation log file
ws.log(msg)
defgo_full():
""" Run the simulation in one go with h.run() """
initialize()
h.run()
# Run
t0 = datetime.now()
# log('Running the simulation in NEURON')
now = datetime.now()
ws.log("Starting the simulation in NEURON on: %s"%now.strftime("%d %h %Y at %H:%M:%S"))
go_by_steps()
now = datetime.now()
ws.log("NEURON simulation finished on: %s"%now.strftime("%d %h %Y at %H:%M:%S"))
log('Whole simulation time: %i seconds'%(now - t0).total_seconds())