// File for generating cell numbers and position values
load_file("Network_distributions.hoc")
load_file("Parameters.hoc")
// Object for position of cells along the -'ve diagonal
objref PR_position
objref PR_population_upper_position[10000], PR_population_lower_position[10000], PR_population_upperright_position[1000]
objref PR_count_upper, PR_count_lower
objref HZ_position
objref HZ_population_upper_position[10000], HZ_population_lower_position[10000], HZ_population_upperright_position[1000]
objref HZ_count_upper, HZ_count_lower
objref BIP_position
objref BIP_population_upper_position[10000], BIP_population_lower_position[10000], BIP_population_upperright_position[1000]
objref BIP_count_upper, BIP_count_lower
objref AII_position
objref AII_population_upper_position[10000], AII_population_lower_position[10000], AII_population_upperright_position[1000]
objref AII_count_upper, AII_count_lower
// Objects for randomly selecting types of cells
objref random_PR, random_BIP
// Calculate cell density values at the distribution
// Arguments are x-postiion, y-position, z-position and cell type index (identifier)
proc calc_cell_density() {
x = $1
y = $2
z = $3
cell_type_index = $4
strdef cell_type
dist_from_fovea = sqrt((x - Fovea_x)^2 + (y - Fovea_y)^2)
if(cell_type_index == 1) {
cell_type = "pr"
cell_density_value = pr_c1*exp(pr_lambda1*dist_from_fovea) + pr_c2*exp(pr_lambda2*dist_from_fovea) + pr_c3*exp(pr_lambda3*dist_from_fovea) + pr_c4*exp(pr_lambda4*dist_from_fovea) + pr_c5*exp(pr_lambda5*dist_from_fovea) + pr_c6*exp(pr_lambda6*dist_from_fovea)
}
if(cell_type_index == 2) {
cell_type = "horizontal"
cell_density_value = horizontal_c1*exp(horizontal_lambda1*dist_from_fovea) + horizontal_c2*exp(horizontal_lambda2*dist_from_fovea) + horizontal_c3*exp(horizontal_lambda3*dist_from_fovea)
}
if(cell_type_index == 3) {
cell_type = "bipolar"
cell_density_value = bipolar_c1*exp(bipolar_lambda1*dist_from_fovea) + bipolar_c2*exp(bipolar_lambda2*dist_from_fovea) + bipolar_c3*exp(bipolar_lambda3*dist_from_fovea) + bipolar_c4*exp(bipolar_lambda4*dist_from_fovea) + bipolar_c5*exp(bipolar_lambda5*dist_from_fovea) + bipolar_c6*exp(bipolar_lambda6*dist_from_fovea) + bipolar_c7*exp(bipolar_lambda7*dist_from_fovea) + bipolar_c8*exp(bipolar_lambda8*dist_from_fovea) + bipolar_c9*exp(bipolar_lambda9*dist_from_fovea)
}
if(cell_type_index == 4) {
cell_type = "AII_amacrine"
cell_density_value = AII_amacrine_c1*exp(AII_amacrine_lambda1*dist_from_fovea) + AII_amacrine_c2*exp(AII_amacrine_lambda2*dist_from_fovea) + AII_amacrine_c3*exp(AII_amacrine_lambda3*dist_from_fovea)
}
if(cell_density_value < 0) {
cell_density_value = 0
}
// print cell_density_value
}
// Calculate the spacing between cells in hex-mosaic based on density values
// Assuming a HCP sphere packing organisation
proc calc_cell_spacing() {
// there are 3 cells per unit lattice, 6 corners containing 1/3 cell and one in the centre
// spacing for a unit lattice
unit_lattice = 3*sqrt(3)/2 // mm^2
stacks = 1
x_shift = 0.5 // mm
y_shift = sqrt(3)/2 // mm
//check whether density = 0 (usually around fovea)
//to account for division by 0, we make the density very small instead
if(cell_density_value < 10) {
cell_density_value = 10
}
cell_density_per_layer = cell_density_value/stacks
cell_spacing_x = x_shift/(sqrt(cell_density_per_layer))
cell_spacing_y = y_shift/(sqrt(cell_density_per_layer))
// print cell_spacing_x
// print cell_spacing_y
}
// Procedure for filling in matrices with positions of cells
proc determine_cell_init_positions() {
PR_position = new Matrix(1000, 2)
HZ_position = new Matrix(1000, 2)
BIP_position = new Matrix(1000, 2)
AII_position = new Matrix(1000, 2)
// -----------Generate Positions for PHOTORECEPTORS------------------------------------------
// generate initial cell at (Eccentricity_x - patch_x/2, Eccentricity_y - patch_y) for PRs (bottom right corner)
calc_cell_density(Eccentricity_x - Patch_size_x/2, Eccentricity_y - Patch_size_y/2, 0, 1)
PR_position.x[0][0] = Eccentricity_x - Patch_size_x/2
PR_position.x[0][1] = Eccentricity_y - Patch_size_y/2
print PR_position.x[0][0]
// generate the diagonal line of cells for populating the patch size
i = 1
while(1) {
calc_cell_density(PR_position.x[i - 1][0], PR_position.x[i - 1][1], 0, 1) // cell_type_index = 1 for PRs
calc_cell_spacing()
PR_position.x[i][0] = PR_position.x[i - 1][0] + cell_spacing_x
PR_position.x[i][1] = PR_position.x[i - 1][1] + cell_spacing_y
if((abs(PR_position.x[i][0] - Eccentricity_x)) > (Patch_size_x/2)) {
PR_total_diag_cells = i + 1
break
}
i = i + 1
}
// -----------Generate Positions for HORIZONTAL CELLS------------------------------------------
// generate initial cell at (Eccentricity_x - patch_x/2, Eccentricity_y - patch_y) for HZs (bottom right corner)
calc_cell_density(Eccentricity_x - Patch_size_x/2, Eccentricity_y - Patch_size_y/2, 0, 2)
HZ_position.x[0][0] = Eccentricity_x - Patch_size_x/2
HZ_position.x[0][1] = Eccentricity_y - Patch_size_y/2
print HZ_position.x[0][0]
// generate the diagonal line of cells for populating the patch size
i = 1
while(1) {
calc_cell_density(HZ_position.x[i - 1][0], HZ_position.x[i - 1][1], 0, 2) // cell_type_index = 2 for HZs
calc_cell_spacing()
HZ_position.x[i][0] = HZ_position.x[i - 1][0] + cell_spacing_x
HZ_position.x[i][1] = HZ_position.x[i - 1][1] + cell_spacing_y
if((abs(HZ_position.x[i][0] - Eccentricity_x)) > (Patch_size_x/2)) {
HZ_total_diag_cells = i + 1
break
}
i = i + 1
}
// -----------Generate Positions for BIPOLAR CELLS------------------------------------------
// generate initial cell at (Eccentricity_x - patch_x/2, Eccentricity_y - patch_y) for BIPs (bottom right corner)
calc_cell_density(Eccentricity_x - Patch_size_x/2, Eccentricity_y - Patch_size_y/2, 0, 3)
BIP_position.x[0][0] = Eccentricity_x - Patch_size_x/2
BIP_position.x[0][1] = Eccentricity_y - Patch_size_y/2
print BIP_position.x[0][0]
// generate the diagonal line of cells for populating the patch size
i = 1
while(1) {
calc_cell_density(BIP_position.x[i - 1][0], BIP_position.x[i - 1][1], 0, 3) // cell_type_index = 3 for BIPs
calc_cell_spacing()
BIP_position.x[i][0] = BIP_position.x[i - 1][0] + cell_spacing_x
BIP_position.x[i][1] = BIP_position.x[i - 1][1] + cell_spacing_y
if((abs(BIP_position.x[i][0] - Eccentricity_x)) > (Patch_size_x/2)) {
BIP_total_diag_cells = i + 1
break
}
i = i + 1
}
// -----------Generate Positions for AII AMACRINE CELLS------------------------------------------
// generate initial cell at (Eccentricity_x - patch_x/2, Eccentricity_y - patch_y) for AIIs (bottom right corner)
calc_cell_density(Eccentricity_x - Patch_size_x/2, Eccentricity_y - Patch_size_y/2, 0, 4)
AII_position.x[0][0] = Eccentricity_x - Patch_size_x/2
AII_position.x[0][1] = Eccentricity_y - Patch_size_y/2
print AII_position.x[0][0]
// generate the diagonal line of cells for populating the patch size
i = 1
while(1) {
calc_cell_density(AII_position.x[i - 1][0], AII_position.x[i - 1][1], 0, 4) // cell_type_index = 4 for AIIs
calc_cell_spacing()
AII_position.x[i][0] = AII_position.x[i - 1][0] + cell_spacing_x
AII_position.x[i][1] = AII_position.x[i - 1][1] + cell_spacing_y
if((abs(AII_position.x[i][0] - Eccentricity_x)) > (Patch_size_x/2)) {
AII_total_diag_cells = i + 1
break
}
i = i + 1
}
}
proc generate_population_positions() {
// -----------Generate Population Positions for PHOTORECEPTORS---------------------------------
PR_count_upper_cells = 0
PR_count_lower_cells = 0
PR_count_upper = new Vector(PR_total_diag_cells)
PR_count_lower = new Vector(PR_total_diag_cells)
// Create positions of upper diagonal portion
for i = 0, PR_total_diag_cells - 1 {
j = 0
PR_population_upper_position[i] = new Matrix(10000, 2)
PR_population_upper_position[i].x[j][0] = PR_position.x[i][0]
PR_population_upper_position[i].x[j][1] = PR_position.x[i][1]
j = 1
while(1) {
calc_cell_density(PR_population_upper_position[i].x[j - 1][0], PR_population_upper_position[i].x[j - 1][1], 0, 1)
calc_cell_spacing()
PR_population_upper_position[i].x[j][0] = PR_population_upper_position[i].x[j - 1][0] - cell_spacing_x
PR_population_upper_position[i].x[j][1] = PR_population_upper_position[i].x[j - 1][1] + cell_spacing_y
if(abs(PR_population_upper_position[i].x[j][0] - Eccentricity_x) > (Patch_size_x/2)) {
PR_count_upper.x[i] = j
PR_count_upper_cells = PR_count_upper_cells + 1
break
}
if(abs(PR_population_upper_position[i].x[j][1] - Eccentricity_y) > (Patch_size_y/2)) {
PR_count_upper.x[i] = j
PR_count_upper_cells = PR_count_upper_cells + 1
break
}
j = j + 1
PR_count_upper_cells = PR_count_upper_cells + 1
}
}
// Create positions of lower diagonal portion
for i = 0, PR_total_diag_cells - 1 {
j = 0
PR_population_lower_position[i] = new Matrix(10000, 2)
PR_population_lower_position[i].x[j][0] = PR_position.x[i][0]
PR_population_lower_position[i].x[j][1] = PR_position.x[i][1]
j = 1
while(1) {
calc_cell_density(PR_population_lower_position[i].x[j - 1][0], PR_population_lower_position[i].x[j - 1][1], 0, 1)
calc_cell_spacing()
PR_population_lower_position[i].x[j][0] = PR_population_lower_position[i].x[j - 1][0] + cell_spacing_x
PR_population_lower_position[i].x[j][1] = PR_population_lower_position[i].x[j - 1][1] - cell_spacing_y
if(abs(PR_population_lower_position[i].x[j][0] - Eccentricity_x) > (Patch_size_x/2)) {
PR_count_lower.x[i] = j
PR_count_lower_cells = PR_count_lower_cells + 1
break
}
j = j + 1
PR_count_lower_cells = PR_count_lower_cells + 1
}
}
// Create final array of positions
// -----------Generate Population Positions for HORIZONTAL CELLS---------------------------------
HZ_count_upper_cells = 0
HZ_count_lower_cells = 0
HZ_count_upperright_cells = 0
HZ_count_upper = new Vector(HZ_total_diag_cells)
HZ_count_lower = new Vector(HZ_total_diag_cells)
// Create positions of upper diagonal portion
for i = 0, HZ_total_diag_cells - 1 {
j = 0
HZ_population_upper_position[i] = new Matrix(10000, 2)
HZ_population_upper_position[i].x[j][0] = HZ_position.x[i][0]
HZ_population_upper_position[i].x[j][1] = HZ_position.x[i][1]
j = 1
while(1) {
calc_cell_density(HZ_population_upper_position[i].x[j - 1][0], HZ_population_upper_position[i].x[j - 1][1], 0, 2)
calc_cell_spacing()
HZ_population_upper_position[i].x[j][0] = HZ_population_upper_position[i].x[j - 1][0] - cell_spacing_x
HZ_population_upper_position[i].x[j][1] = HZ_population_upper_position[i].x[j - 1][1] + cell_spacing_y
if(abs(HZ_population_upper_position[i].x[j][0] - Eccentricity_x) > (Patch_size_x/2)) {
HZ_count_upper.x[i] = j
HZ_count_upper_cells = HZ_count_upper_cells + 1
break
}
if(abs(HZ_population_upper_position[i].x[j][1] - Eccentricity_y) > (Patch_size_y/2)) {
HZ_count_upper.x[i] = j
HZ_count_upper_cells = HZ_count_upper_cells + 1
break
}
j = j + 1
HZ_count_upper_cells = HZ_count_upper_cells + 1
}
}
// Create positions of lower diagonal portion
for i = 0, HZ_total_diag_cells - 1 {
j = 0
HZ_population_lower_position[i] = new Matrix(10000, 2)
HZ_population_lower_position[i].x[j][0] = HZ_position.x[i][0]
HZ_population_lower_position[i].x[j][1] = HZ_position.x[i][1]
j = 1
while(1) {
calc_cell_density(HZ_population_lower_position[i].x[j - 1][0], HZ_population_lower_position[i].x[j - 1][1], 0, 2)
calc_cell_spacing()
HZ_population_lower_position[i].x[j][0] = HZ_population_lower_position[i].x[j - 1][0] + cell_spacing_x
HZ_population_lower_position[i].x[j][1] = HZ_population_lower_position[i].x[j - 1][1] - cell_spacing_y
if(abs(HZ_population_lower_position[i].x[j][0] - Eccentricity_x) > (Patch_size_x/2)) {
HZ_count_lower.x[i] = j
HZ_count_lower_cells = HZ_count_lower_cells + 1
break
}
j = j + 1
HZ_count_lower_cells = HZ_count_lower_cells + 1
}
}
// -----------Generate Population Positions for BIPOLAR CELLS---------------------------------
BIP_count_upper_cells = 0
BIP_count_lower_cells = 0
BIP_count_upperright_cells = 0
BIP_count_upper = new Vector(BIP_total_diag_cells)
BIP_count_lower = new Vector(BIP_total_diag_cells)
// Create positions of upper diagonal portion
for i = 0, BIP_total_diag_cells - 1 {
j = 0
BIP_population_upper_position[i] = new Matrix(10000, 2)
BIP_population_upper_position[i].x[j][0] = BIP_position.x[i][0]
BIP_population_upper_position[i].x[j][1] = BIP_position.x[i][1]
j = 1
while(1) {
calc_cell_density(BIP_population_upper_position[i].x[j - 1][0], BIP_population_upper_position[i].x[j - 1][1], 0, 3)
calc_cell_spacing()
BIP_population_upper_position[i].x[j][0] = BIP_population_upper_position[i].x[j - 1][0] - cell_spacing_x
BIP_population_upper_position[i].x[j][1] = BIP_population_upper_position[i].x[j - 1][1] + cell_spacing_y
if(abs(BIP_population_upper_position[i].x[j][0] - Eccentricity_x) > (Patch_size_x/2)) {
BIP_count_upper.x[i] = j
BIP_count_upper_cells = BIP_count_upper_cells + 1
break
}
if(abs(BIP_population_upper_position[i].x[j][1] - Eccentricity_y) > (Patch_size_y/2)) {
BIP_count_upper.x[i] = j
BIP_count_upper_cells = BIP_count_upper_cells + 1
break
}
j = j + 1
BIP_count_upper_cells = BIP_count_upper_cells + 1
}
}
// Create positions of lower diagonal portion
for i = 0, BIP_total_diag_cells - 1 {
j = 0
BIP_population_lower_position[i] = new Matrix(10000, 2)
BIP_population_lower_position[i].x[j][0] = BIP_position.x[i][0]
BIP_population_lower_position[i].x[j][1] = BIP_position.x[i][1]
j = 1
while(1) {
calc_cell_density(BIP_population_lower_position[i].x[j - 1][0], BIP_population_lower_position[i].x[j - 1][1], 0, 3)
calc_cell_spacing()
BIP_population_lower_position[i].x[j][0] = BIP_population_lower_position[i].x[j - 1][0] + cell_spacing_x
BIP_population_lower_position[i].x[j][1] = BIP_population_lower_position[i].x[j - 1][1] - cell_spacing_y
if(abs(BIP_population_lower_position[i].x[j][0] - Eccentricity_x) > (Patch_size_x/2)) {
BIP_count_lower.x[i] = j
BIP_count_lower_cells = BIP_count_lower_cells + 1
break
}
j = j + 1
BIP_count_lower_cells = BIP_count_lower_cells + 1
}
}
// -----------Generate Population Positions for AII AMACRINE CELLS---------------------------------
AII_count_upper_cells = 0
AII_count_lower_cells = 0
AII_count_upperright_cells = 0
AII_count_upper = new Vector(AII_total_diag_cells)
AII_count_lower = new Vector(AII_total_diag_cells)
// Create positions of upper diagonal portion
for i = 0, AII_total_diag_cells - 1 {
j = 0
AII_population_upper_position[i] = new Matrix(10000, 2)
AII_population_upper_position[i].x[j][0] = AII_position.x[i][0]
AII_population_upper_position[i].x[j][1] = AII_position.x[i][1]
j = 1
while(1) {
calc_cell_density(AII_population_upper_position[i].x[j - 1][0], AII_population_upper_position[i].x[j - 1][1], 0, 4)
calc_cell_spacing()
AII_population_upper_position[i].x[j][0] = AII_population_upper_position[i].x[j - 1][0] - cell_spacing_x
AII_population_upper_position[i].x[j][1] = AII_population_upper_position[i].x[j - 1][1] + cell_spacing_y
if(abs(AII_population_upper_position[i].x[j][0] - Eccentricity_x) > (Patch_size_x/2)) {
AII_count_upper.x[i] = j
AII_count_upper_cells = AII_count_upper_cells + 1
break
}
if(abs(AII_population_upper_position[i].x[j][1] - Eccentricity_y) > (Patch_size_y/2)) {
AII_count_upper.x[i] = j
AII_count_upper_cells = AII_count_upper_cells + 1
break
}
j = j + 1
AII_count_upper_cells = AII_count_upper_cells + 1
}
}
// Create positions of lower diagonal portion
for i = 0, AII_total_diag_cells - 1 {
j = 0
AII_population_lower_position[i] = new Matrix(10000, 2)
AII_population_lower_position[i].x[j][0] = AII_position.x[i][0]
AII_population_lower_position[i].x[j][1] = AII_position.x[i][1]
j = 1
while(1) {
calc_cell_density(AII_population_lower_position[i].x[j - 1][0], AII_population_lower_position[i].x[j - 1][1], 0, 4)
calc_cell_spacing()
AII_population_lower_position[i].x[j][0] = AII_population_lower_position[i].x[j - 1][0] + cell_spacing_x
AII_population_lower_position[i].x[j][1] = AII_population_lower_position[i].x[j - 1][1] - cell_spacing_y
if(abs(AII_population_lower_position[i].x[j][0] - Eccentricity_x) > (Patch_size_x/2)) {
AII_count_lower.x[i] = j
AII_count_lower_cells = AII_count_lower_cells + 1
break
}
j = j + 1
AII_count_lower_cells = AII_count_lower_cells + 1
}
}
}
func calc_PR_density() {
x = $1
y = $2
dist_from_fovea = sqrt((x - Fovea_x)^2 + (y - Fovea_y)^2)
rod_density = rod_c1*exp(rod_lambda1*dist_from_fovea) + rod_c2*exp(rod_lambda2*dist_from_fovea) + rod_c3*exp(rod_lambda3*dist_from_fovea)
cone_density = cone_c1*exp(cone_lambda1*dist_from_fovea) + cone_c2*exp(cone_lambda2*dist_from_fovea) + cone_c3*exp(cone_lambda3*dist_from_fovea)
total_PR_density = int(rod_density + cone_density)
random_PR = new Random(total_PR_density)
d = random_PR.uniform(0, total_PR_density)
if(d < rod_density) {
PR_type_value = 1 // 1 for rod
} else {
PR_type_value = 2 // 2 for cone
}
return PR_type_value
}
func calc_BIP_density() {
x = $1
y = $2
dist_from_fovea = sqrt((x - Fovea_x)^2 + (y - Fovea_y)^2)
rod_bipolar_density = rod_bipolar_c1*exp(rod_bipolar_lambda1*dist_from_fovea) + rod_bipolar_c2*exp(rod_bipolar_lambda2*dist_from_fovea) + rod_bipolar_c3*exp(rod_bipolar_lambda3*dist_from_fovea)
ON_bipolar_density = ON_bipolar_c1*exp(ON_bipolar_lambda1*dist_from_fovea) + ON_bipolar_c2*exp(ON_bipolar_lambda2*dist_from_fovea) + ON_bipolar_c3*exp(ON_bipolar_lambda3*dist_from_fovea)
OFF_bipolar_density = OFF_bipolar_c1*exp(OFF_bipolar_lambda1*dist_from_fovea) + OFF_bipolar_c2*exp(OFF_bipolar_lambda2*dist_from_fovea) + OFF_bipolar_c3*exp(OFF_bipolar_lambda3*dist_from_fovea)
total_BIP_density = int(rod_bipolar_density + ON_bipolar_density + OFF_bipolar_density)
random_BIP = new Random(total_BIP_density)
for f = 0, Seed {
d = random_BIP.uniform(0, total_BIP_density)
}
if(d < rod_bipolar_density) {
BIP_type_value = 1 // 1 for rod bipolar
}
if(d >= rod_bipolar_density && d < (rod_bipolar_density + ON_bipolar_density)) {
BIP_type_value = 2 // 2 for ON bipolar
}
if(d >= (rod_bipolar_density + ON_bipolar_density)) {
BIP_type_value = 3 // 3 for OFF bipolar
}
return BIP_type_value
}