Skip to content
This repository has been archived by the owner on Jul 18, 2024. It is now read-only.

Axon class corrected #4

Closed
wants to merge 11 commits into from
Closed

Axon class corrected #4

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
3a99c0b
unstage my additional directories
arashgmn May 15, 2021
d6f1315
some refactoring
arashgmn May 27, 2021
a09bd83
Update Launcher_OSS_lite.py
arashgmn Jun 17, 2021
0af3482
Update MRI_DTI_prep_new.py
arashgmn Jun 17, 2021
2ba15a9
Merge pull request #1 from SFB-ELAINE/master
arashgmn Jun 17, 2021
b8d2c66
corrected the Axon class for McIntyre's model
arashgmn Jun 23, 2021
5722f40
Merge branch 'SFB-ELAINE:master' into master
arashgmn Jun 23, 2021
0fd0591
more refactoring
arashgmn Jun 23, 2021
2a8ea43
update gitignore to skip pycache folders
arashgmn Jun 23, 2021
39be161
restoring to the default input_dict
arashgmn Jun 23, 2021
9336d9f
fixed assymetric pattern in Neural_array_processing + updated gitignore
arashgmn Jun 24, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions .gitignore
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,2 +1,4 @@
access-token.txt
local-build.sh
OSS_platform/**/__pycache__

Binary file modified OSS-DBS - tutorial.pdf
View file
Open in desktop
Binary file not shown.
86 changes: 53 additions & 33 deletions OSS_platform/Axon_files/axon.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -138,21 +138,21 @@ def get_axonparams(self):
'total_nodes' int
Total number of compartments
'ranvier_nodes' int
Number of node of Ranvier compartments
Number of node of Ranvier compartments
'para1_nodes' int
Number of first paranodal compartments
Number of first paranodal compartments (MYSA)
'para2_nodes' int
Number of second paranodal compartments
Number of second paranodal compartments (FLUT)
'inter_nodes' int
Number of internodal compartments
Number of internodal compartments (STIN)
'ranvier_length' float
Length of the Node of Ranvier
'para1_length' float
Length of the first paranodal compartments
Length of the first paranodal compartments (MYSA)
'para2_length' float
Length of the second paranodal compartments
Length of the second paranodal compartments (FLUT)
'inter_length' float
Length of the internodal compartments
Length of the internodal compartments (STIN)
'deltax' float
Length from a node of Ranvier to the next
'fiberD' float
Expand All @@ -162,9 +162,9 @@ def get_axonparams(self):
'node_diameter': float
Diameter of the nodes of Ranvier
'para1_diameter': float
Diameter of the first paranodal compartments
Diameter of the first paranodal compartments (MYSA)
'para2_diameter': float,
Diameter of the second paranodal compartments
Diameter of the second paranodal compartments (FLUT)
'condg': float
Axial conductivity
'lamellas': int
Expand All @@ -191,23 +191,40 @@ def __create_ref_nodes(self):
ncomp = nstin+nmysa+nflut+1

ref_nodes = np.zeros(n,)
for i in range(0, n):
j = i % ncomp
if j == 0 or j == 1: # mysa <-> ranvier node
if i == 0: # first of all nodes
ref_nodes[i] = l_ranvier/2.
else:
ref_nodes[i] = ref_nodes[i-1]+l_para1/2.+l_ranvier/2.
elif (j > 1 and j < (int(nmysa/2)+1)) or (j > (int(nmysa/2)+nflut+nstin) and j <(nmysa+nflut+nstin)): # mysa <-> mysa node
ref_nodes[i] = ref_nodes[i-1]+l_para1
elif j == (int(nmysa/2)+1) or j == (int(nmysa/2)+nflut+nstin): # mysa <-> flut node
ref_nodes[i] = ref_nodes[i-1]+l_para1/2.+l_para2/2.
elif (j > (int(nmysa/2)+1) and j < (int(nmysa/2)+int(nflut/2)+1)) or (j > (int(nmysa/2)+int(nflut/2)+nstin) and j < (int(nmysa/2)+nflut+nstin)): # flut <-> flut node
ref_nodes[i] = ref_nodes[i-1]+l_para2
elif j == (int(nmysa/2)+int(nflut/2)+1) or j == (int(nmysa/2)+int(nflut/2)+nstin): # flut <-> stin node
ref_nodes[i] = ref_nodes[i-1]+l_para2/2.+l_inter/2.
else: # stin <-> stin node
ref_nodes[i] = ref_nodes[i-1]+l_inter

# An easier implementation based on the structure of the internodals segments
structure = 'RMF'+6*'S'+'FM' # Node of Rnavier + internodal structure

# Let's rename the lengths, just for convenience
L_R = l_ranvier
L_F = l_para2
L_M = l_para1
L_S = l_inter

compartment = [0] # intial node of Ranvier
for i in range(n-1):
l1 = eval('L_'+structure[i%len(structure)]) # current segment length
l2 = eval('L_'+structure[(i+1)%len(structure)]) # next segment length
compartment.append(compartment[-1]+(l1+l2)/2)

ref_nodes = np.array(compartment)
# for i in range(0, n):
# j = i % ncomp
# if j == 0 or j == 1: # mysa <-> ranvier node
# if i == 0: # first of all nodes
# ref_nodes[i] = l_ranvier/2.
# else:
# ref_nodes[i] = ref_nodes[i-1]+l_para1/2.+l_ranvier/2.
# elif (j > 1 and j < (int(nmysa/2)+1)) or (j > (int(nmysa/2)+nflut+nstin) and j <(nmysa+nflut+nstin)): # mysa <-> mysa node
# ref_nodes[i] = ref_nodes[i-1]+l_para1
# elif j == (int(nmysa/2)+1) or j == (int(nmysa/2)+nflut+nstin): # mysa <-> flut node
# ref_nodes[i] = ref_nodes[i-1]+l_para1/2.+l_para2/2.
# elif (j > (int(nmysa/2)+1) and j < (int(nmysa/2)+int(nflut/2)+1)) or (j > (int(nmysa/2)+int(nflut/2)+nstin) and j < (int(nmysa/2)+nflut+nstin)): # flut <-> flut node
# ref_nodes[i] = ref_nodes[i-1]+l_para2
# elif j == (int(nmysa/2)+int(nflut/2)+1) or j == (int(nmysa/2)+int(nflut/2)+nstin): # flut <-> stin node
# ref_nodes[i] = ref_nodes[i-1]+l_para2/2.+l_inter/2.
# else: # stin <-> stin node
# ref_nodes[i] = ref_nodes[i-1]+l_inter

return ref_nodes * 1e-6 # um -> m

Expand All @@ -217,6 +234,9 @@ def __create_nodes(self):
self.__ref_nodes[-1]/2
else:
ref_nodes = self.__ref_nodes

# The reference nodes are not 100% symmetric, i.e. sometimes ref_nodes[i]!=ref_nodes[-i]
# due to the fitting option
# standard position(along x axis in [x,y,z])
nodes = np.zeros((len(ref_nodes),3))
nodes[:,0] = ref_nodes
Expand Down Expand Up @@ -256,17 +276,17 @@ def __create_axon_parameters(self, diameter):

axon_parameters = copy.deepcopy(self.AXONPARAMETERS[diameter])

nranvier = axon_parameters['ranvier_nodes']
nstin = int(float(axon_parameters['inter_nodes'])/(nranvier-1))
nmysa = int(float(axon_parameters['para1_nodes'])/(nranvier-1))
nflut = int(float(axon_parameters['para2_nodes'])/(nranvier-1))
nranvier = axon_parameters['ranvier_nodes'] #21
nstin = int(float(axon_parameters['inter_nodes'])/(nranvier-1)) # 6
nmysa = int(float(axon_parameters['para1_nodes'])/(nranvier-1)) # 2
nflut = int(float(axon_parameters['para2_nodes'])/(nranvier-1)) # 2

axon_parameters['fiberD']=diameter
axon_parameters['total_nodes'] = axon_parameters['ranvier_nodes'] + \
axon_parameters['para1_nodes']+axon_parameters['para2_nodes'] + \
axon_parameters['inter_nodes']
axon_parameters['inter_nodes'] #221
axon_parameters['inter_length'] = (axon_parameters['deltax'] -
axon_parameters['ranvier_length'] -
nmysa*axon_parameters['para1_length'] -
nflut*axon_parameters['para2_length'])/float(nstin)
return axon_parameters
nflut*axon_parameters['para2_length'])/float(nstin) #70.5
return axon_parameters
7 changes: 4 additions & 3 deletions OSS_platform/CSF_refinement_new.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -304,6 +304,7 @@ def Refine_CSF(MRI_param,DTI_param,Scaling,Domains,Field_calc_param,rel_div,CSF_
z_neuron_min=min(Vertices_neur[:,2])

#go over all voxels and check whether it contains CSF and intersect with the mesh
eps = 0.000000001
for x_coord in x_vect:
for y_coord in y_vect:
for z_coord in z_vect:
Expand All @@ -314,9 +315,9 @@ def Refine_CSF(MRI_param,DTI_param,Scaling,Domains,Field_calc_param,rel_div,CSF_

if (x_pos<=x_neuron_max+CSF_ref_add and x_pos>=x_neuron_min-CSF_ref_add and y_pos<=y_neuron_max+CSF_ref_add and y_pos>=y_neuron_min-CSF_ref_add and z_pos<=z_neuron_max+CSF_ref_add and z_pos>=z_neuron_min-CSF_ref_add):

xv_mri=int((x_coord)/MRI_param.x_vox_size-0.000000001) #defines number of steps to get to the voxels containing x[0] coordinate
yv_mri=(int((y_coord)/MRI_param.y_vox_size-0.000000001))*MRI_param.M_x #defines number of steps to get to the voxels containing x[0] and x[1] coordinates
zv_mri=(int((z_coord)/MRI_param.z_vox_size-0.000000001))*MRI_param.M_x*MRI_param.M_y #defines number of steps to get to the voxels containing x[0], x[1] and x[2] coordinates
xv_mri= x_coord//(MRI_param.x_vox_size-eps) #defines number of steps to get to the voxels containing x[0] coordinate
yv_mri=(y_coord//(MRI_param.y_vox_size-eps))*MRI_param.M_x #defines number of steps to get to the voxels containing x[0] and x[1] coordinates
zv_mri=(z_coord//(MRI_param.z_vox_size-eps))*MRI_param.M_x*MRI_param.M_y #defines number of steps to get to the voxels containing x[0], x[1] and x[2] coordinates
glob_index=int(xv_mri+yv_mri+zv_mri)

pnt=Point(x_pos,y_pos,z_pos)
Expand Down
71 changes: 54 additions & 17 deletions OSS_platform/FEM_in_spectrum.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -198,8 +198,12 @@ class Conductivity : public dolfin::Expression

if int(sine_freq)==int(signal_freq):
c_unscaled = CompiledExpression(compile_cpp_code(conductivity_code).Conductivity(),
c00=c_unscaled00, c01=c_unscaled01, c02=c_unscaled02, c11=c_unscaled11, c12=c_unscaled12, c22=c_unscaled22, degree=0)
tensor= as_tensor([[c_unscaled[0], c_unscaled[1], c_unscaled[2]], [c_unscaled[1], c_unscaled[3], c_unscaled[4]],[c_unscaled[2],c_unscaled[4],c_unscaled[5]]])
c00=c_unscaled00, c01=c_unscaled01, c02=c_unscaled02,
c11=c_unscaled11, c12=c_unscaled12, c22=c_unscaled22,
degree=0)
tensor = as_tensor([[c_unscaled[0], c_unscaled[1], c_unscaled[2]],
[c_unscaled[1], c_unscaled[3], c_unscaled[4]],
[c_unscaled[2],c_unscaled[4],c_unscaled[5]]])
f_vector_repr=project(tensor,TensorFunctionSpace(mesh, "Lagrange", 1),solver_type="cg", preconditioner_type="amg")
file=File('Tensors/Ellipsoids_unscaled_at_'+str(signal_freq)+'_Hz.pvd')
file<<f_vector_repr
Expand Down Expand Up @@ -423,8 +427,16 @@ def solve_Laplace(Sim_setup,Solver_type,Vertices_array,Domains,core,Full_IFFT,ou
Sim_setup.mesh.coordinates()
Sim_setup.mesh.init()

# to get conductivity (and permittivity if EQS formulation) mapped accrodingly to the subdomains. k_val_r is just a list of conductivities (S/mm!) in a specific order to scale the cond. tensor
kappa,k_val_r=get_dielectric_properties_from_subdomains(Sim_setup.mesh,Sim_setup.subdomains,Sim_setup.Laplace_eq,Domains.Float_contacts,Sim_setup.conductivities,Sim_setup.rel_permittivities,Sim_setup.sine_freq)
# to get conductivity (and permittivity if EQS formulation) mapped accrodingly to the subdomains.
# k_val_r is just a list of conductivities (S/mm!) in a specific order to scale the cond. tensor
kappa,k_val_r = get_dielectric_properties_from_subdomains(Sim_setup.mesh,
Sim_setup.subdomains,
Sim_setup.Laplace_eq,
Domains.Float_contacts,
Sim_setup.conductivities,
Sim_setup.rel_permittivities,
Sim_setup.sine_freq)

if int(Sim_setup.sine_freq)==int(Sim_setup.signal_freq):
file=File('Field_solutions/Conductivity_map_'+str(Sim_setup.signal_freq)+'Hz.pvd')
file<<kappa[0]
Expand All @@ -434,19 +446,33 @@ def solve_Laplace(Sim_setup,Solver_type,Vertices_array,Domains,core,Full_IFFT,ou

# to get tensor scaled by the conductivity map
if Sim_setup.anisotropy==1:
Cond_tensor=get_scaled_cond_tensor(Sim_setup.mesh,Sim_setup.subdomains,Sim_setup.sine_freq,Sim_setup.signal_freq,Sim_setup.unscaled_tensor,k_val_r)
Cond_tensor = get_scaled_cond_tensor(Sim_setup.mesh,
Sim_setup.subdomains,
Sim_setup.sine_freq,
Sim_setup.signal_freq,
Sim_setup.unscaled_tensor,
k_val_r)
else:
Cond_tensor=False #just to initialize

#In case of current-controlled stimulation, Dirichlet_bc or the whole potential distribution will be scaled afterwards (due to the system's linearity)
V_space,Dirichlet_bc,ground_index=get_solution_space_and_Dirichlet_BC(Sim_setup.c_c,Sim_setup.mesh,Sim_setup.boundaries,Sim_setup.element_order,Sim_setup.Laplace_eq,Domains.Contacts,Domains.fi)
Cond_tensor = False #just to initialize

# In case of current-controlled stimulation, Dirichlet_bc
# the whole potential distribution will be scaled afterwards
# (due to the system's linearity)
V_space, Dirichlet_bc, ground_index = get_solution_space_and_Dirichlet_BC(Sim_setup.c_c,
Sim_setup.mesh,
Sim_setup.boundaries,
Sim_setup.element_order,
Sim_setup.Laplace_eq,
Domains.Contacts,
Domains.fi)
#ground index refers to the ground in .med/.msh file

# to solve the Laplace equation div(kappa*grad(phi))=0 (variational form: a(u,v)=L(v))
phi_sol=define_variational_form_and_solve(V_space,Dirichlet_bc,kappa,Sim_setup.Laplace_eq,Cond_tensor,Solver_type)
phi_sol =define_variational_form_and_solve(V_space, Dirichlet_bc, kappa, Sim_setup.Laplace_eq,
Cond_tensor, Solver_type)

if Sim_setup.Laplace_eq=='EQS':
(phi_r,phi_i)=phi_sol.split(deepcopy=True)
(phi_r,phi_i)= phi_sol.split(deepcopy=True)
else:
phi_r=phi_sol
phi_i=Function(V_space)
Expand All @@ -460,7 +486,8 @@ def solve_Laplace(Sim_setup,Solver_type,Vertices_array,Domains,core,Full_IFFT,ou


if Sim_setup.c_c==1 or Sim_setup.CPE_status==1: #we compute E-field, currents and impedances only for current-controlled or if CPE is used
J_ground=get_current(Sim_setup.mesh,Sim_setup.boundaries,Sim_setup.element_order,Sim_setup.Laplace_eq,Domains.Contacts,kappa,Cond_tensor,phi_r,phi_i,ground_index)
J_ground = get_current(Sim_setup.mesh, Sim_setup.boundaries, Sim_setup.element_order, Sim_setup.Laplace_eq,
Domains.Contacts, kappa, Cond_tensor, phi_r, phi_i, ground_index)
#If EQS, J_ground is a complex number

#V_across=max(Domains.fi[:], key=abs) # voltage drop in the system
Expand All @@ -484,25 +511,34 @@ def solve_Laplace(Sim_setup,Solver_type,Vertices_array,Domains,core,Full_IFFT,ou
print("Currently, CPE can be used only for simulations with two contacts. Please, assign the rest to 'None'")
raise SystemExit

Dirichlet_bc_with_CPE,total_impedance=get_CPE_corrected_Dirichlet_BC(Sim_setup.boundaries,Sim_setup.CPE_param,Sim_setup.Laplace_eq,Sim_setup.sine_freq,Sim_setup.signal_freq,Domains.Contacts,Domains.fi,V_across,Z_tissue,V_space)
Dirichlet_bc_with_CPE, total_impedance = get_CPE_corrected_Dirichlet_BC(Sim_setup.boundaries,
Sim_setup.CPE_param,
Sim_setup.Laplace_eq,
Sim_setup.sine_freq,
Sim_setup.signal_freq,
Domains.Contacts, Domains.fi,
V_across, Z_tissue, V_space)

f=open('Field_solutions/Impedance'+str(core)+'.csv','ab')
np.savetxt(f, total_impedance, delimiter=" ")
f.close()

# to solve the Laplace equation for the adjusted Dirichlet
phi_sol_CPE=define_variational_form_and_solve(V_space,Dirichlet_bc_with_CPE,kappa,Sim_setup.Laplace_eq,Cond_tensor,Solver_type)
phi_sol_CPE = define_variational_form_and_solve(V_space, Dirichlet_bc_with_CPE, kappa,
Sim_setup.Laplace_eq, Cond_tensor, Solver_type)
if Sim_setup.Laplace_eq=='EQS':
(phi_r_CPE,phi_i_CPE)=phi_sol_CPE.split(deepcopy=True)
else:
phi_r_CPE=phi_sol_CPE
phi_i_CPE=Function(V_space)
phi_i_CPE.vector()[:] = 0.0

J_ground_CPE=get_current(Sim_setup.mesh,Sim_setup.boundaries,Sim_setup.element_order,Sim_setup.Laplace_eq,Domains.Contacts,kappa,Cond_tensor,phi_r_CPE,phi_i_CPE,ground_index)
J_ground_CPE = get_current(Sim_setup.mesh, Sim_setup.boundaries, Sim_setup.element_order,
Sim_setup.Laplace_eq, Domains.Contacts, kappa, Cond_tensor,
phi_r_CPE, phi_i_CPE, ground_index)

# just resaving
phi_sol,phi_r,phi_i,J_ground=(phi_sol_CPE,phi_r_CPE,phi_i_CPE,J_ground_CPE)
phi_sol, phi_r, phi_i, J_ground = (phi_sol_CPE, phi_r_CPE, phi_i_CPE, J_ground_CPE)

if Full_IFFT==1:
Hdf=HDF5File(Sim_setup.mesh.mpi_comm(), "Field_solutions_functions/solution"+str(np.round(Sim_setup.sine_freq,6))+".h5", "w")
Expand Down Expand Up @@ -567,7 +603,8 @@ def solve_Laplace(Sim_setup,Solver_type,Vertices_array,Domains,core,Full_IFFT,ou
else:
Dirichlet_bc_scaled.append(DirichletBC(V_space, np.real((Domains.fi[bc_i])/J_ground),Sim_setup.boundaries,Domains.Contacts[bc_i]))

phi_sol_check=define_variational_form_and_solve(V_space,Dirichlet_bc_scaled,kappa,Sim_setup.Laplace_eq,Cond_tensor,Solver_type)
phi_sol_check = define_variational_form_and_solve(V_space, Dirichlet_bc_scaled, kappa,
Sim_setup.Laplace_eq, Cond_tensor, Solver_type)

if Sim_setup.Laplace_eq=='EQS':
(phi_r_check,phi_i_check)=phi_sol_check.split(deepcopy=True)
Expand Down
Loading