%% This code is written by Nooshin Abdollahi - Oct 23 2022
% Matlab codes necessary for peripheral network modeling.
% The output of this code will be saved and loaded in Python script.
clc
clear all
close all
%%
R = [4;4] ; % two fibers with exactly the same radius.
C(1,1) = 0 ;
C(1,2) = 0 ;
C(2,1) = 8.1 ; % the edge to edge distance between the two fibers is 4 microns. (12-4*2) = 4
C(2,2) = 0 ;
G(1,1) = 4 ; % Both fibers are Abeta fibers (group 4) ;
G(1,2) = 4 ;
figure
for i=1:length(R)
viscircles(C(i,:),R(i,1), 'color', 'k','Linewidth',0.2)
hold on
end
% xlim([-20 20])
% ylim([-20 20])
axis equal
%%
unique_radius = unique(R) ;
%% To find morphological parameters
i = 1 ;
x= unique_radius(i,1)*2 ; % fiber diameter
p1 = 0.01876 ;
p2 = 0.4787 ;
p3 = 0.1204 ;
axonD = p1*x^2 + p2*x + p3 ;
paraD2 = axonD ;
p1 = 0.006304 ;
p2 = 0.2071 ;
p3 = 0.5339 ;
nodeD = p1*x^2 + p2*x + p3 ;
paraD1 = nodeD ;
p1 = -8.215 ;
p2 = 272.4 ;
p3 = -780.2 ;
deltax = p1*x^2 + p2*x + p3 ;
p1 = -0.0199 ;
p2 = 3.016 ;
p3 = 17.44 ;
paralength2 = p1*x^2 + p2*x + p3 ; % FLUT length
paralength2_segments = 5 ;
p1 = -0.389 ;
p2 = 14.88 ;
p3 = 9.721 ;
nl = round ( p1*x^2 + p2*x + p3 ) ;
paralength1=3 ;
nodelength=1.0 ;
interlength=(deltax-nodelength-(2*paralength1)-(2*paralength2)) ;
internode_segments = 40 ; % number of segments in one internode
parameters(i,1) = axonD ;
parameters(i,2) = nodeD ;
parameters(i,3) = deltax ;
parameters(i,4) = paralength2/paralength2_segments ;
parameters(i,5) = nl ;
parameters(i,6) = interlength/internode_segments ; % length of each segment
%%
N=zeros(length(C),length(C)) ;
dist=zeros(length(C),length(C)) ;
Distance_edge=zeros(length(C),length(C)) ;
Distance_Centers = pdist2(C, C) ; % distance between centers
%%%%%%% Distance between edges of circles not centers of circles %%%%%%%%%
for i=1:length(Distance_Centers)
for j=1:length(Distance_Centers)
if i ~= j
Distance_edge(i,j)= Distance_Centers(i,j) - (R(i,1)+R(j,1)) ;
end
end
end
%%
m= 1000 ; % minimum distance (threshold)
for k=1:length(Distance_Centers)
for i=1:length(Distance_Centers)
if k ~= i
if Distance_Centers(k,i)< m
N(k,i)=1 ; % 1 means connection, 0 means no connection
dist(k,i) = Distance_Centers(k,i) ; % distance between the axons
end
end
end
end
for j=1:length(N)
for i=1:length(N)
if N(j,i)==1
N(i,j)=0 ;
dist(i,j) = 0 ;
end
end
end
%%
Total_Length_Axon = 31500 ; % the total length of axons is "about" ... (microns)
%%%%%%%%%%%%%%%%%%%%% the other type of axons %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%% there is only one type
for J=1: length(unique_radius)
z1=0 ;
z2=0 ;
z3=0 ;
z4=0 ;
j=1 ;
i=2 ;
which_radius = J ;
%
eval(['coordinates_sections_' num2str(J) '(1,1) = 0 ;']);
eval(['coordinates_sections_' num2str(J) '(1,i) = coordinates_sections_' num2str(J) '(1,i-1)+nodelength;']);
eval(['mid_coordinates_sections_' num2str(J) '(1,j) = 0.5;']);
eval(['name_sections_' num2str(J) '(1,j) = sprintf("node_%d",z1) ;']);
z1=z1+1 ;
i=i+1 ;
j=j+1 ;
eval(['coordinates_sections_' num2str(J) '(1,i) = coordinates_sections_' num2str(J) '(1,i-1)+paralength1;']);
eval(['mid_coordinates_sections_' num2str(J) '(1,j) = (coordinates_sections_' num2str(J) '(1,i) + coordinates_sections_' num2str(J) '(1,i-1) )/2;']);
eval(['name_sections_' num2str(J) '(1,j) = sprintf("MYSA_%d",z2) ;']);
z2=z2+1 ;
i=i+1 ;
j=j+1 ;
for rr=1:paralength2_segments
eval(['coordinates_sections_' num2str(J) '(1,i) = coordinates_sections_' num2str(J) '(1,i-1)+parameters(which_radius,4);']);
eval(['mid_coordinates_sections_' num2str(J) '(1,j) = (coordinates_sections_' num2str(J) '(1,i) + coordinates_sections_' num2str(J) '(1,i-1) )/2;']);
eval(['name_sections_' num2str(J) '(1,j) = sprintf("FLUT_%d",z3) ;']);
z3=z3+1 ;
i=i+1 ;
j=j+1 ;
end
for r=1:internode_segments
eval(['coordinates_sections_' num2str(J) '(1,i) = coordinates_sections_' num2str(J) '(1,i-1)+(parameters(which_radius,6));']);
eval(['mid_coordinates_sections_' num2str(J) '(1,j) = (coordinates_sections_' num2str(J) '(1,i) + coordinates_sections_' num2str(J) '(1,i-1) )/2;']);
eval(['name_sections_' num2str(J) '(1,j) = sprintf("STIN_%d",z4) ;']);
z4=z4+1 ;
i=i+1 ;
j=j+1 ;
end
for rr=1:paralength2_segments
eval(['coordinates_sections_' num2str(J) '(1,i) = coordinates_sections_' num2str(J) '(1,i-1)+parameters(which_radius,4);']);
eval(['mid_coordinates_sections_' num2str(J) '(1,j) = (coordinates_sections_' num2str(J) '(1,i) + coordinates_sections_' num2str(J) '(1,i-1) )/2;']);
eval(['name_sections_' num2str(J) '(1,j) = sprintf("FLUT_%d",z3) ;']);
z3=z3+1 ;
i=i+1 ;
j=j+1 ;
end
eval(['coordinates_sections_' num2str(J) '(1,i) = coordinates_sections_' num2str(J) '(1,i-1)+paralength1;']);
eval(['mid_coordinates_sections_' num2str(J) '(1,j) = (coordinates_sections_' num2str(J) '(1,i) + coordinates_sections_' num2str(J) '(1,i-1) )/2;']);
eval(['name_sections_' num2str(J) '(1,j) = sprintf("MYSA_%d",z2) ;']);
z2=z2+1 ;
i=i+1 ;
j=j+1 ;
Length_one_set(J,1) = eval(['coordinates_sections_' num2str(J) '(1,i-1)']) ; % Length of one set including (node+MYSA+5*FLUT+40*Internodes+5*FLUT+MYSA)
number_sets(J,1) = round (Total_Length_Axon/Length_one_set(J,1)) ; % how many sets we need for this specific radius in order to get a total length about 10000 microns
%
for ii=2:number_sets(J,1) % it starts from 2 becasue the first set is already above
eval(['coordinates_sections_' num2str(J) '(1,i) = coordinates_sections_' num2str(J) '(1,i-1)+nodelength;']);
eval(['mid_coordinates_sections_' num2str(J) '(1,j) = (coordinates_sections_' num2str(J) '(1,i) + coordinates_sections_' num2str(J) '(1,i-1) )/2;']);
eval(['name_sections_' num2str(J) '(1,j) = sprintf("node_%d",z1) ;']);
z1=z1+1 ;
i=i+1 ;
j=j+1 ;
eval(['coordinates_sections_' num2str(J) '(1,i) = coordinates_sections_' num2str(J) '(1,i-1)+paralength1;']);
eval(['mid_coordinates_sections_' num2str(J) '(1,j) = (coordinates_sections_' num2str(J) '(1,i) + coordinates_sections_' num2str(J) '(1,i-1) )/2;']);
eval(['name_sections_' num2str(J) '(1,j) = sprintf("MYSA_%d",z2) ;']);
z2=z2+1 ;
i=i+1 ;
j=j+1 ;
for rr=1:paralength2_segments
eval(['coordinates_sections_' num2str(J) '(1,i) = coordinates_sections_' num2str(J) '(1,i-1)+parameters(which_radius,4);']);
eval(['mid_coordinates_sections_' num2str(J) '(1,j) = (coordinates_sections_' num2str(J) '(1,i) + coordinates_sections_' num2str(J) '(1,i-1) )/2;']);
eval(['name_sections_' num2str(J) '(1,j) = sprintf("FLUT_%d",z3) ;']);
z3=z3+1 ;
i=i+1 ;
j=j+1 ;
end
for r=1:internode_segments
eval(['coordinates_sections_' num2str(J) '(1,i) = coordinates_sections_' num2str(J) '(1,i-1)+(parameters(which_radius,6));']);
eval(['mid_coordinates_sections_' num2str(J) '(1,j) = (coordinates_sections_' num2str(J) '(1,i) + coordinates_sections_' num2str(J) '(1,i-1) )/2;']);
eval(['name_sections_' num2str(J) '(1,j) = sprintf("STIN_%d",z4) ;']);
z4=z4+1 ;
i=i+1 ;
j=j+1 ;
end
for rr=1:paralength2_segments
eval(['coordinates_sections_' num2str(J) '(1,i) = coordinates_sections_' num2str(J) '(1,i-1)+parameters(which_radius,4);']);
eval(['mid_coordinates_sections_' num2str(J) '(1,j) = (coordinates_sections_' num2str(J) '(1,i) + coordinates_sections_' num2str(J) '(1,i-1) )/2;']);
eval(['name_sections_' num2str(J) '(1,j) = sprintf("FLUT_%d",z3) ;']);
z3=z3+1 ;
i=i+1 ;
j=j+1 ;
end
eval(['coordinates_sections_' num2str(J) '(1,i) = coordinates_sections_' num2str(J) '(1,i-1)+paralength1;']);
eval(['mid_coordinates_sections_' num2str(J) '(1,j) = (coordinates_sections_' num2str(J) '(1,i) + coordinates_sections_' num2str(J) '(1,i-1) )/2;']);
eval(['name_sections_' num2str(J) '(1,j) = sprintf("MYSA_%d",z2) ;']);
z2=z2+1 ;
i=i+1 ;
j=j+1 ;
end
Length_All_sets(J,1) = eval(['coordinates_sections_' num2str(J) '(1,i-1)']) ; % "around" 10000 microns
Number_of_nodes (J,1) = z1 ;
end
%%
figure
plot(coordinates_sections_1 ,1 , 'ob')
hold on
plot(coordinates_sections_1 ,1.05 , 'or')
ylim([1 1.5])
%%
compartments_in_a_set = 3 + internode_segments + (2*paralength2_segments) ;
%%
%%%%%%%%%%%%%%%%%%%%% Connect each axon to itself %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for s=1 :length((unique_radius))
for t=1 :length((unique_radius))
if t== s
clear S_first
clear S_second
clear n nn f1 f2 fz m mm
L1=0 ;
L2=0 ;
p=1 ;
for n=0:Number_of_nodes(s,1)-1
S_first(p,1) = sprintf("node_%d",n) ;
S_second(p,1) = sprintf("node_%d",n) ;
p=p+1 ;
end
eval(['Connect_types_' num2str(s),num2str(t) '(1:length(S_first),1) = S_first ;']);
eval(['Connect_types_' num2str(s),num2str(t) '(length(S_first)+1:length(S_first)+length(S_second),1) = S_second ;']);
%%%%%%%%%% save data %%%%%%%%%%
clear SAVE
SAVE = char( eval(['Connect_types_' num2str(s),num2str(t)])) ;
filename=['Connect_types_' num2str(s),num2str(t) '.mat'];
save(filename , 'SAVE');
end
end
end
%% Output
% *** Connection_types_xx are saved above
save('R.txt' , 'R' , '-ascii') % save radius of all axons
save('C.txt' , 'C' , '-ascii') % save location of all axons (x,y)
save('G.txt' , 'G' , '-ascii') % save Groups of all axons. i.e., each radius in "R" belongs to which type (Group1,2,3 or 4)
% To be used in NEURON to build models
save('parameters.txt' , 'parameters' , '-ascii')
save('Number_of_nodes.txt' , 'Number_of_nodes' , '-ascii')
save('unique_radius.txt' , 'unique_radius' , '-ascii')
save('neighboringAxon.txt' , 'N' , '-ascii')
% % save('neighbouringAxon.mat','N')
save('dist.txt' , 'dist' , '-ascii')
save('Distance_edge.txt' , 'Distance_edge' , '-ascii')
%% Next: Boundary Condition
% finding the neurons on the edge, for now I do it manually
neurons_on_edge_up = [0 1 ] ;
% neurons_on_edge_right = [ ] ;
% neurons_on_edge_down = [ ] ;
% neurons_on_edge_left = [ ] ;
neurons_on_edge = [neurons_on_edge_up] ;
% finding the boundary cable coordinates
distance_boundary_to_neuron = 7 ; % microns
for i=1:length(neurons_on_edge_up)
B_coordinates_up(i,1) = C(neurons_on_edge_up(i)+1,1) ;
B_coordinates_up(i,2) = C(neurons_on_edge_up(i)+1,2) + distance_boundary_to_neuron ;
end
Boundary_coordinates(:,1) = [B_coordinates_up(:,1)'] ;
Boundary_coordinates(:,2) = [B_coordinates_up(:,2)'] ;
Boundary_coordinates ;
plot(Boundary_coordinates(:,1),Boundary_coordinates(:,2),'o','Color',[0.9290 0.6940 0.1250])
%% finding the neighboring neurons to which boundary cables are connected
C_and_Boundary_coordinates(:,1) = [ C(:,1)' , Boundary_coordinates(:,1)' ] ;
C_and_Boundary_coordinates(:,2) = [ C(:,2)' , Boundary_coordinates(:,2)' ] ;
size_C = size(C) ;
size_Boundary_coordinates = size(Boundary_coordinates) ;
Boundary_N=zeros(length(C_and_Boundary_coordinates),length(C_and_Boundary_coordinates)) ;
Boundary_dist=zeros(length(C_and_Boundary_coordinates),length(C_and_Boundary_coordinates)) ;
Boundary_Distance = pdist2(C_and_Boundary_coordinates, C_and_Boundary_coordinates) ;
m= 8 ; % minimum distance (threshold)
for k=1:length(Boundary_Distance)
for i=1:length(Boundary_Distance)
if k ~= i
if Boundary_Distance(k,i)< m
Boundary_N(k,i)=1 ;
Boundary_dist(k,i)=Boundary_Distance(k,i) ;
end
end
end
end
Boundary_neighboring = Boundary_N(size_C+1:size(Boundary_N) ,1:size_C) ;
Boundary_dist = Boundary_dist(size_C+1:size(Boundary_N) ,1:size_C) ;
%% find the sections for connections
% for Boundary Cable
z=0 ;
i=2 ;
j=1 ;
number_of_sections = 4000 ;
sectionlength= Total_Length_Axon/number_of_sections ;
for ii=0:number_of_sections-1
Boundary_coordinates_sections(1,1) = 0 ;
Boundary_coordinates_sections(1,i) = Boundary_coordinates_sections(1,i-1) + sectionlength ;
Boundary_mid_coordinates_section(1,j) = (Boundary_coordinates_sections(1,i) + Boundary_coordinates_sections(1,i-1) )/2 ;
Boundary_name_sections(1,j) = sprintf("section_%d",z) ;
i=i+1 ;
z=z+1 ;
j=j+1 ;
end
for t=1 :length((unique_radius))
clear S_first
clear S_second
clear n f1 m
L=0 ;
p=1 ;
for n=1:compartments_in_a_set:compartments_in_a_set*Number_of_nodes(t,1)
f1 = eval(['mid_coordinates_sections_' num2str(t) '(1,n)']) ;
for m=1:length (Boundary_coordinates_sections)-1
if f1 >= Boundary_coordinates_sections(1,m) && f1 <= Boundary_coordinates_sections(1,m+1)
S_first(p,1) = sprintf("node_%d",L) ;
S_second(p,1) = Boundary_name_sections(1,m) ;
p=p+1 ;
end
end
L=L+1 ;
end
eval(['Boundary_to_' num2str(t) '(1:length(S_first),1) = S_first ;']);
eval(['Boundary_to_' num2str(t) '(length(S_first)+1:length(S_first)+length(S_second),1) = S_second ;']);
%%%%%%%%%% save data %%%%%%%%%%
clear SAVE
SAVE = char( eval(['Boundary_to_' num2str(t)])) ;
filename=['Boundary_to_' num2str(t) '.mat'];
save(filename , 'SAVE');
end
%%% Boundary to Boundary are not connected
%% Output
% *** Boundary_to_x are saved above
save('Boundary_coordinates.txt' , 'Boundary_coordinates' , '-ascii') % (x,y)
save('Boundary_neighboring.txt' , 'Boundary_neighboring' , '-ascii')
save('Boundary_dist.txt' , 'Boundary_dist' , '-ascii')
%% Next step: Transverse resistances (g) unit: S/cm2
% Calculating c in PLOS computational biology paper, Equation 9
for i=1 :length((unique_radius))
for j=1 :length((unique_radius))
if j>=i
clear Rg
z=2 ;
L = length(eval(['Connect_types_' num2str(i),num2str(j)]))/2 ;
Rg=zeros(L,1) ;
for Z=2:L-1
CO1 = find( eval(['name_sections_' num2str(i)]) == eval(['Connect_types_' num2str(i),num2str(j) '(Z-1,1)']) );
CO2 = find( eval(['name_sections_' num2str(i)]) == eval(['Connect_types_' num2str(i),num2str(j) '(Z+1,1)']));
g1 = eval(['mid_coordinates_sections_' num2str(i) '(1,CO1)']) ;
g2 = eval(['mid_coordinates_sections_' num2str(i) '(1,CO2)']) ;
Rg(z,1) = (g2-g1)/2 ;
z=z+1 ;
if g2 < g1
msg = 'Error occurred. g2<g1';
error(msg)
end
end
CO1 = find( eval(['name_sections_' num2str(i)]) == eval(['Connect_types_' num2str(i),num2str(j) '(1,1)']) );
CO2 = find( eval(['name_sections_' num2str(i)]) == eval(['Connect_types_' num2str(i),num2str(j) '(2,1)']));
g1 = eval(['mid_coordinates_sections_' num2str(i) '(1,CO1)']) ;
g2 = eval(['mid_coordinates_sections_' num2str(i) '(1,CO2)']) ;
Rg(1,1) = (g2-g1)/2 ;
CO1 = find( eval(['name_sections_' num2str(i)]) == eval(['Connect_types_' num2str(i),num2str(j) '(L-1,1)']) );
CO2 = find( eval(['name_sections_' num2str(i)]) == eval(['Connect_types_' num2str(i),num2str(j) '(L,1)']));
g1 = eval(['mid_coordinates_sections_' num2str(i) '(1,CO1)']) ;
g2 = eval(['mid_coordinates_sections_' num2str(i) '(1,CO2)']) ;
Rg(L,1) = (g2-g1)/2 ;
%%%%%%%%%% save data %%%%%%%%%%
filename=['Rg_' num2str(i),num2str(j) '.txt'];
save(filename , 'Rg' , '-ascii')
end
end
end
%%
%%% Areas of the nodes
for i=1 :length((unique_radius))
for j=2 :length((unique_radius))
if j>i
clear Areas
k=0 ;
L = length(eval(['Connect_types_' num2str(i),num2str(j)]))/2 ;
Areas =zeros(L,1) ;
for Z=1:L
if (eval(['Connect_types_' num2str(i),num2str(j) '(Z,1)']) == sprintf("node_%d",k) )
Areas(Z,1) = pi*parameters(i+1,2) ; % surface area = pi*d*L (L=1), parameters(i,2): nodeD
k=k+1 ;
else
Areas(Z,1) = pi*parameters(j+1,2) ; % surface area = pi*d*L (L=1), parameters(i,2): nodeD
end
end
%%%%%%%%%% save data %%%%%%%%%%
filename=['Areas_' num2str(i),num2str(j) '.txt'];
save(filename , 'Areas' , '-ascii')
end
end
end
%% between boundary cables and the axons %%%%%
for j=1:length((unique_radius))
clear Rg
z=2 ;
L = length(eval(['Boundary_to_' num2str(j)]))/2 ;
Rg=zeros(L,1) ;
for Z=2:L-1
CO1 = find( eval(['name_sections_' num2str(j)]) == eval(['Boundary_to_' num2str(j) '(Z-1,1)']) );
CO2 = find( eval(['name_sections_' num2str(j)]) == eval(['Boundary_to_' num2str(j) '(Z+1,1)']) );
g1 = eval(['mid_coordinates_sections_' num2str(j) '(1,CO1)']) ;
g2 = eval(['mid_coordinates_sections_' num2str(j) '(1,CO2)']) ;
Rg(z,1) = (g2-g1)/2 ;
z=z+1 ;
end
CO1 = find( eval(['name_sections_' num2str(j)]) == eval(['Boundary_to_' num2str(j) '(1,1)']) );
CO2 = find( eval(['name_sections_' num2str(j)]) == eval(['Boundary_to_' num2str(j) '(2,1)']));
g1 = eval(['mid_coordinates_sections_' num2str(j) '(1,CO1)']) ;
g2 = eval(['mid_coordinates_sections_' num2str(j) '(1,CO2)']) ;
Rg(1,1) = (g2-g1)/2 ;
CO1 = find( eval(['name_sections_' num2str(j)]) == eval(['Boundary_to_' num2str(j) '(L-1,1)']) );
CO2 = find( eval(['name_sections_' num2str(j)]) == eval(['Boundary_to_' num2str(j) '(L,1)']));
g1 = eval(['mid_coordinates_sections_' num2str(j) '(1,CO1)']) ;
g2 = eval(['mid_coordinates_sections_' num2str(j) '(1,CO2)']) ;
Rg(L,1) = (g2-g1)/2 ;
%%%%%%%%%% save data %%%%%%%%%%
filename=['Boundary_Rg_' num2str(j) '.txt'];
save(filename , 'Rg' , '-ascii')
end
%%
b=1 ;
for Z=1:(length(Connect_types_11)/2)-1
COO1 = find( eval(['name_sections_' num2str(1)]) == eval(['Connect_types_' num2str(1),num2str(1) '(Z+1,1)']) );
COO2 = find( eval(['name_sections_' num2str(1)]) == eval(['Connect_types_' num2str(1),num2str(1) '(Z,1)']));
gg1 = eval(['mid_coordinates_sections_' num2str(i) '(1,COO1)']) ;
gg2 = eval(['mid_coordinates_sections_' num2str(i) '(1,COO2)']) ;
nodes_dist(b,1) = gg1-gg2 ;
b=b+1 ;
end