# Python program using Lennard-Jones force to approximate cortical stem-cell close packing symmetry

The program repeatedly goes from a squeeze cycle to align points into symmetric as possible geometry, and expansion cycle:

I’m still new to Python, so please forgive where it may be obvious that I need to figure out how to most simply improve array structure and other things. Coding suggestions are appreciated:

``````# Form stable vesicle membrane around sphere,  Gary Gaulin, 8/4/2019
# Uses Lennard Jones Potential for attraction and repulsion behavior,
# representative of polar molecules, bubbles, and cortical cell migration.
# Starts with almost enough points to fill sphere.
# There are then large holes and unfilled 5 sided pores.
# For cortical stem cells these pores are a space for them to fill, not be any.
# Radius is then slowly decreased until pressure causes chain reaction of realignments,
# towards most stable configuration, ideal symmetry, for given number of points.

import sys
import pygame
import random
import numpy as np

Pi = 3.14159265358979323846264338327950288412
HalfPi = Pi / 2

# Startup Parameters
# Example of 392 is the dimple count of a vintage Hogan golfball.
# It has an unusually uniform dimple uniformity, as though made by packing real balls around it.
# Other examples are:
# points = 12     # https://en.wikipedia.org/wiki/Regular_icosahedron
# points = 32     # https://en.wikipedia.org/wiki/Pentakis_dodecahedron
points = 392                        # Total number of (cells, molecules, etc.) points.
screensize = Width, Height = (640, 640)
sas = 4*Pi*(sr*sr)                  # Area of sphere at this given screen radius.
sar = sas/points                    # Divide for area of one Lennard-Jones radius, between points.
r = np.sqrt(sar/Pi)*2               # Radius between points, given area.
spacing = r-((1/points)*r)          # Radius between points, ideal spacing.
spheremaxradius = int(sr * 1.2)     # Maximum radius of sphere, points are located around.
pressuremax = 20                    # Maximum amount of average Lennard-Jones force to apply.
pressurestep = 1                    # Amount of radius change per time step, in screen pixels.
pressuretime = 200                  # Time step duration to hold under pressure before releasing.
dt = .01                            # Timestep size, how (slower but) precise it moves.
maxvelocity = spacing/10            # Velocity limit for when enough force to fly off sphere.
cutoffdist = spacing*3              # Distance two points are too far away to exert force.
connectdist = spacing*1.15          # Distance two points are touching, bumping each other.
cutoffcalctime = 50                 # Time steps between updates of i,j points within cut-off range.
rotationstep = .05                  # Amount of rotation per timestep.
rotation = 0                        # Initialize rotation angle, to draw all points on screen.
timestepnumber = 0

epsilon = 1                         # Lennard Jones variable, remains 1, normalized.
sigma = 1                           # Lennard Jones variable, remains 1, normalized.

Blue = (0, 0, 150)
Black = (0, 0, 0)
Green = (0, 100, 0)
Yellow = (120, 120, 0)

# Initialize x,y,z arrays for points on sphere.
xp = [0 for _ in range(points)]
yp = [0 for _ in range(points)]
zp = [0 for _ in range(points)]
# Initialize Lennard-Jones sum array.
ljsum = [0 for _ in range(points)]
# X,Y of points on screen.
xyzs = [0 for _ in range(points)]
# Periodically updated list of points within cutoff range.
near = [0 for _ in range(points*100)]
# Periodically updated list of connections between points.
connect = [0 for _ in range(points*20)]

# Cartesian to Spherical coordinates.
def cart2sphere(x, y, z):
hxy = np.hypot(x, y)
el = np.arctan2(z, hxy)
az = np.arctan2(y, x)
return az, el, np.hypot(hxy, z)

# Spherical to Cartesian coordinates.
def sphere2cart(az, el, ra):
theta = ra * np.cos(el)
x = theta * np.cos(az)
y = theta * np.sin(az)
z = ra * np.sin(el)
return x, y, z

def drawaxisline(x, y, z, c):
az, el, r = cart2sphere(x, y, z)  # Cartesian to Spherical (az, el, r)
x1, y1, z1 = sphere2cart(az + rotation, el, r)  # Spherical to Cartesian.
pygame.draw.line(screen, c, (xc, yc),  (x1+xc, z1+yc), 1)

def findnearby():
ni = 0
ci = 0
for i in range(points):             # For all the points in array
for j in range(points):         # Except itself
if i != j:
dx = xp[j] - xp[i]      # Calculate x,y,z distances.
dy = yp[j] - yp[i]
dz = zp[j] - zp[i]
# Use the three distances for calculating an always positive 3D distance.
# This is the very time consuming command to less occasionally use.
dist = np.sqrt((dx * dx) + (dy * dy) + (dz * dz))
# If distance is within cutoff range then save 'j' point, at this 'i'.
if dist < cutoffdist:
near[ni] = i, j
ni += 1
# If distance is within touching range then save 'j' point, at this 'i'.
if dist < connectdist:
connect[ci] = i, j
ci += 1
print(ni, ci)
return ni, ci

def update_positions():
# Initialize velocity vector arrays for summing forces, combined angles and magnitudes.
xi = [0 for _ in range(points)]
yi = [0 for _ in range(points)]
zi = [0 for _ in range(points)]
xj = [0 for _ in range(points)]
yj = [0 for _ in range(points)]
zj = [0 for _ in range(points)]
# Clear sum , and maximum
ljt = 0                 # Initialize variable for totaling Lennard-Jones force.
ljmax = -10000          # Initialize maximum variable, negative when points attract.
# For all the points in array.
for n in range(ni):
i,j = near[n]
# Calculate x,y,z distances.
dx = xp[j] - xp[i]
dy = yp[j] - yp[i]
dz = zp[j] - zp[i]
# Use the three distances for calculating an always positive 3D distance.
dist = np.sqrt((dx * dx) + (dy * dy) + (dz * dz))
# Calculate Lennard Jones force.
ljr = dist/spacing              # Lennard-Jones radius. At ljr=1 force is 0.
if ljr > 0:                     # Prevent divide by zero error when LJ radius is zero.
lj = 4 * epsilon * ((sigma/ljr)**12 - (sigma/ljr)**6)
else:
lj = 0
# Update monitor variables for totaling average and maximum Lennard-Jones reading..
if lj > ljmax:
ljmax = lj
ljt += lj                       # Add to sum for average of all points.
ljsum[i] += lj                  # Used for brightness to draw on screen.
# Reduce timestep for smooth movement, fewer points wildly bouncing off each other.
lj = lj * dt
# Limit amount of displacement/velocity as though through water, or cell migration, not frictionless vacuum.
if lj > maxvelocity:
lj = maxvelocity
if lj < -maxvelocity:
lj = -maxvelocity
# Use distance for calculating spherical azimuth and elevation angle r, as in cart2sphere.
el = np.arctan2(dz, dist)
az = np.arctan2(dy, dx)
# Use Lennard Jones amount to calculate displacement amount, for given force.
x, y, z = sphere2cart(az, el, lj)
# Update x,y,z force acting on i point, from this j point.
xi[i] -= x
yi[i] -= y
zi[i] -= z
# Update other point, j. For every action there is equal and opposite reaction.
xj[j] += x
yj[j] += y
zj[j] += z
# Add force displacement amounts, to points.
for i in range(points):
xp[i] += xi[i] + xj[i]
yp[i] += yi[i] + yj[i]
zp[i] += zi[i] + zj[i]
return ljt/ni, ljmax

def render():
screen.fill(Black)
# Precalculate all the rotated X,Y screen locations and their colors.
for i in range(points):
az, el, r = cart2sphere(xp[i], yp[i], zp[i])    # Cartesian to Spherical (az, el, r)
x, y, z = sphere2cart(az + rotation, el, r)     # Spherical to Cartesian.
# Vary brightness by amount of force, positive increasingly orange, negative blue.
poscolor = 0
negcolor = 0
if ljsum[i] > 0:
poscolor = int(ljsum[i])
if poscolor > 255:
poscolor = 255
if ljsum[i] < 0:
negcolor = int(-ljsum[i] * 50)
if negcolor > 255:
negcolor = 255
c = (poscolor, (poscolor + negcolor) * .6, negcolor)
xyzs[i] = x, y, z, c
# To properly draw on top of each other sort the list according to what is furthest away, [1] = by Y axis.
xyzs.sort(key=lambda x: x[1])
# Draw all points to screen.
for i in range(points):
# Calculate x,y and radius then draw point to screen buffer.
x, y, z, c = xyzs[i]
xs = int(x + xc)
ys = int(z + yc)
d = ((y+sphereradius)/spherewidth*.4) + .6   # Fractional amount, to 1, according to how far away, depth.
if d > 1:
d = 1
cr = int(d * spacing * .45)             # Circle radius.
pygame.draw.circle(screen, c, [xs, ys], cr, 0)
pygame.draw.circle(screen, (255, 255, 255), [xs, ys], cr, 1)
# Draw axis lines.
ra = sphereradius * 1.1                 # Distance of axis line, from 0.
drawaxisline(0, 0, ra, Yellow)
if rotation > ThreeQtrRadian:           # Green x behind blue y?
drawaxisline(ra, 0, 0, Green)
drawaxisline(0, ra, 0, Blue)
if rotation < ThreeQtrRadian:           # Green x on top of blue y?
drawaxisline(ra, 0, 0, Green)
pygame.display.update()                     # Display everything drawn into screen buffer.
fps.tick(60)

# Fill Environment arrays with randomly placed points.
# They will become located on surface of sphere.
for i in range(points):
xp[i] = random.randrange(-100, 100)         # Fill positive half of cube.
yp[i] = random.randrange(-100, 100)
zp[i] = random.randrange(-100, 100)

fnd = cutoffcalctime          # Less than full amount when starting off.
dialog = "Squeezing"
holdcount = pressuretime

# PyGame Initialization
pygame.init()
screen = pygame.display.set_mode(screensize)
xc = int(screensize[0]/2)
yc = int(screensize[1]/2)
pygame.display.set_caption("Spherical Membrane.   Positive Axis: x=green, y=blue, z=yellow")
fps = pygame.time.Clock()
paused = False

while True:
for event in pygame.event.get():
if event.type == pygame.QUIT:
pygame.quit()
sys.exit()
if event.type == pygame.KEYUP:
if event.key == pygame.K_SPACE:
paused = not paused
if not paused:
timestepnumber += 1
# Keep points on sphere by converting to Spherical then (using sphere radius) to Cartesian again.
for i in range(points):  # 1D Radial distance between x points.
az, el, r = cart2sphere(xp[i], yp[i], zp[i])  # Cartesian to Spherical.
xp[i], yp[i], zp[i] = sphere2cart(az, el, sphereradius)  # Spherical to Cartesian again.
# Update positions of points, and rotation to drawn on screen.
fnd += 1
if fnd > cutoffcalctime:
fnd = 0
ni, ci = findnearby()
ljsum = [0 for _ in range(points)]
ljavr, ljmax = update_positions()
render()
rotation += rotationstep                    # Add to screen rotation amount.
if rotation > Radian:                       # Keep in range from 0 to 1 Radian
rotation = 0
# Pressurize and unpressurize by changing radius of sphere, by average L-J force amount.
moreless = ""                               # Clear optional "More" or "Less" dialog.
if dialog == "Squeezing":
holdcount -= 1                          # Decrement squeezing time counter.
if holdcount < 0:                       # Time to recalculate points within cutoff range?
dialog = "Begin Expansion LJ=0"     # Changing 'dialog' ends Squeezing.
if ljavr <pressuremax:
moreless = "More, "                 # Add to dialog as "Squeezing More".
else:
moreless = "Less, "                 # Add to dialog as "Squeezing Less".
else:
if ljavr < 0:                           # Fully expanded?
# Save file here.
dialog = "Squeezing"                # Changing dialog to "Squeezing" ends expansion.
holdcount = pressuretime            # Start new countdown to next expansion cycle.
else:
print(timestepnumber, "  ", dialog, moreless, "Radius =", int(sphereradius), "  LJaverage =", ljavr, "  LJmax =", ljmax)``````

My next challenge is to use this to articulate a network, to attempt propagating traveling waves in the usual simple way I use, around a 3D brain ball or surface with some curvature as in biological membranes.

And credit to my wife Laurie, for counting all the dimples in ten’s using a magic marker, while I worked on the code to recreate the pattern. It’s sometimes hard to force out the influence of a small amount of space, golf ball dimple mould engineers could have in other ways squeezed out too.

5 Likes

Wow!

When it seemed Laurie might say no I told her I would give her credit for helping. To make sure it’s the right number required a couple of recounts including one where I marked off the dots (now with official count of 392 instead of 393) and helped with pictures to show the unusually regular hexagonal pattern. I just edited the code in the opening post to include number change.

Out of a good sized selection my Mom long ago found around the next door golf course it was the only one with impressive hexagonality. After my having tried geodesic dome websites and other possibilities this is still the best test example I could find. Previously I was wasting time guessing what to try next, so being within at least one away from almost perfect symmetry was for me a breakthrough.

The ball you showed is close to the perfection I was looking for, but would have been rejected. All others had different sized dimples, were too shallow to have been impressed by spheres, a mostly square pattern, or too much empty space to be a good example. There was though this one exception, using 1950’s technology, where to me it looks like thousands of man hours may have been spent close packing exactly the right number of beads or bearings to the most ideal/stable geometry possible.

This code should eliminate all need for an earlier mentioned hexagonal math methodology, only good for flat surfaces anyway. It’s skipping an in-between 2D hexagonal axis system based representation that leaves out much detail, by going straight to a 3D x,y,z location in a virtual brain of any shape. This can seem like getting more fancy than necessary but starting with articulation files (I will next have the program generate) makes it all a what connects to where problem where waves of activity already go from region to region then hopefully back to where it originated from, without having to go from one 2D network to another as separate entities.

As in a real cortex (especially mouse and other unfolded types) there is a continuous sheet to subdivide into specialized areas, where at least in visual processing horizontal traveling waves were recorded passing through multiple regions. Question is then what can that activity alone do. For a spherical brain ball it’s in a way a looking through straw type view that might be pieced together from their timing and signal direction. A wave started by itself is expected to arrive back from all directions at the same time, while from others the wave would have a delay in the way all directions of it orderly pass through.

3 Likes

To be confident of the code I ended up needing to improve the physics related math. Regardless of what I did the irregularities remained, but after some time this nice 5 and 6 connection symmetry emerged:

Contacting 6 neighbors is shown in green, 5 in yellow.

``````# View files from Vesicle Maker program, Gary Gaulin, 12/25/2019

import sys
import pygame
import random
import numpy as np

Pi = 3.14159265358979323846264338327950288412
HalfPi = Pi / 2

Blue = (0, 0, 200)
Black = (0, 0, 0)
Green = (0, 100, 0)
Yellow = (250, 250, 0)
Red = (255, 0, 0)

rotationstep = .005                         # Amount of rotation per timestep.
rotation = 0                                # Initialize rotation angle, for drawing all points on screen.
timestepnumber = 0                          # For displaying the time step number.
filesavedcount = 0                          # For numbering the XYZPoints_?.txt files.
filecount = 0

# File to read must have no blank lines or comments added to it, comma separated values only.
destination = []                        # Clear destination list lists are appended to.
with open(filename, 'r') as f:          # Open file in read mode using 'f' as its name:
for line in f:                      # For each comma separated line of code in 'f':
s = line.rstrip('\n')
d = []  # Empty list named 'd' to fill with elements.
if len(s) > 0:
s = s.split(',')             # Split line at commas into list named 's'
for e in s:                     # For each string element listed in 's' list:
d.append(float(e))          # Append floating point element value to list.
destination.append(d)           # Append list of elements to destination list.
return destination              # Returns list of lists only = [[,,],[,,],[,,]]

filename1 = 'xyz' + str(filenumber) + 'Points.txt'
filename2 = 'xyz' + str(filenumber) + 'Connect.txt'
dialog = 'Loaded ' + filename1 + ' and ' + filename2
try:
f = open(filename1)
f.close()
except IOError:
else:
global xyzpoints
global connect
print(dialog)
return dialog

# Cartesian to Spherical coordinates.
def cart2sphere(xyz):
hxy = np.hypot(xyz[0], xyz[1])
el = np.arctan2(xyz[2], hxy)
az = np.arctan2(xyz[1], xyz[0])
return az, el, np.hypot(hxy, xyz[2])

# Spherical to Cartesian coordinates.
def sphere2cart(az, el, ra):
theta = ra * np.cos(el)
x = theta * np.cos(az)
y = theta * np.sin(az)
z = ra * np.sin(el)
return x, y, z

def drawaxisline(x, y, z, c):
az, el, r = cart2sphere([x, y, z])  # Cartesian to Spherical (az, el, r)
x1, y1, z1 = sphere2cart(az + rotation, el, r)  # Spherical to Cartesian.
pygame.draw.line(screen, c, (xc, yc),  (x1+xc, z1+yc), 1)

def update_connections():
maxconnectspace = xyzpoints[pointcount]
for i in range(pointcount):  # For all the i points
for j in range(pointcount):  # and for all the j points
if i != j:  # Except itself
dx = xyzpoints[j][0] - xyzpoints[i][0]  # Calculate x,y,z distances.
dy = xyzpoints[j][1] - xyzpoints[i][1]
dz = xyzpoints[j][2] - xyzpoints[i][2]
# Use the three distances for calculating an always positive 3D distance.
dist = np.sqrt((dx * dx) + (dy * dy) + (dz * dz))
# If distance is within cutoff range then save 'j' point, at this 'i'.
if dist < maxconnectspace:
connect[i].append(j)
return connect

def render():
screen.fill(Black)
# Precalculate all the rotated X,Y screen locations and their colors.
for i in range(pointcount):
az, el, r = cart2sphere(xyzpoints[i])           # Cartesian to Spherical (az, el, r)
x, y, z = sphere2cart(az + rotation, el, r)     # Spherical to Cartesian.
# Vary brightness by amount of force, positive increasingly orange, negative blue.
if len(connect[i]) == 6:
c = Green
else:
if len(connect[i]) == 5:
c = Yellow
else:
c = Red
xyzscreen[i] = x, y, z, c
# To properly draw on top of each other sort the list according to what is furthest away, [1] = by Y axis.
xyzscreen.sort(key=lambda x: x[1])
# Draw all points to screen.
for i in range(pointcount):
# Calculate x,y and radius then draw point to screen buffer.
x, y, z, c = xyzscreen[i]
xs = int(x + xc)
ys = int(z + yc)
cr = int(spacing/2)
pygame.draw.circle(screen, c, [xs, ys], cr, 0)
pygame.draw.circle(screen, (255, 255, 255), [xs, ys], cr, 1)
pygame.display.update()                     # Display everything drawn into screen buffer.
fps.tick(60)

# Start program
filecount += 1
pointcount = len(xyzpoints)
xyzscreen = [0 for _ in range(pointcount)]      # X,Y to draw on screen.
screensize = Width, Height = (640, 640)
sar = sas/pointcount                        # Divide for area of one Lennard-Jones radius, between points.
rbp = np.sqrt(sar/Pi)*2                     # Radius between points, for given area.
spacing = rbp-((1.2/pointcount)*rbp)        # Radius between points, spacing to also fit screen.
print("Spacing", spacing,rbp)

# PyGame Initialization
pygame.init()
screen = pygame.display.set_mode(screensize)
xc = int(screensize[0]/2)
yc = int(screensize[1]/2)
pygame.display.set_caption("Vesicle. Left and right arrow changes file.")
fps = pygame.time.Clock()
paused = False

while True:
for event in pygame.event.get():
if event.type == pygame.QUIT:
pygame.quit()
sys.exit()
if event.type == pygame.KEYUP:
if event.key == pygame.K_SPACE:
paused = not paused
if event.key == pygame.K_RIGHT:
filecount += 1
filecount -= 1
if event.key == pygame.K_LEFT:
filecount -= 1
filecount += 1
if not paused:
timestepnumber += 1
# Print everything to screen. Will include any movement off the sphere during last time step.
render()
rotation += rotationstep                    # Add to screen rotation amount.
if rotation > Radian:                       # Keep in range from 0 to 1 Radian
rotation = 0``````

While being squeezed there are visible hotspots where 1 of 6 gets squeezed out, in order to wrap the flat 2D network geometry around the sphere:

``````# Make stable vesicle membrane around sphere,  Gary Gaulin, 12/25/2019
# Attraction/repulsion force is calculated using Lennard Jones Potential.
# Particles representative of polar molecules, bubbles, and cortical cell migration
# are put into motion using the Velocity Verlet integrator.
#  https://discourse.numenta.org/t/python-program-using-lennard-jones-force-to-approximate-cortical-stem-cell-close-packing-symmetry/6456

import sys
import pygame
import random
import numpy as np

Pi = 3.14159265358979323846264338327950288412
HalfPi = Pi / 2
Blue = (0, 0, 150)
Black = (0, 0, 0)
Green = (0, 120, 0)
Yellow = (100, 100, 0)
Red = (255, 0, 0)
White = (255, 255, 255)

# Startup Parameters
# Test example of 392 is from dimple count of a 1950's Ben Hogan-5 golfball.
# It has an unusually uniform dimple uniformity, as though mould packed 292 real balls together.
pointcount = 392
screensize = Width, Height = (640, 640)     # This is for PyGame window vesicle is shown inside.
screenfillcolor = Black                     # Black background. Later writing of file makes white, to flash when saved.
sar = sas/pointcount                        # Divide for area of one Lennard-Jones radius, between points.
rbp = np.sqrt(sar/Pi)*2                     # Radius between points, for given area.
spacing = rbp-((1.2/pointcount)*rbp)        # Radius between points, spacing to also fit screen.
print("Spacing", spacing, rbp)              # First display ideal spacing in pixels between points.
connectdist = spacing*1.2                   # Distance two points are touching, bumping each other.
cutoffdist = spacing*2.5                    # Distance two points are too far away to influence each other.
cutoffcalctime = 40                         # Time steps between updates of i,j points within cut-off range.
rotationstep = .05                          # Amount of rotation per timestep.
rotation = 0                                # Initialize rotation angle, for drawing all points on screen.
pressurestep = spacing/50                   # Amount of radius change per time step, in screen pixels.
pressuremax = 80                            # Maximum amount of average Lennard-Jones force to apply.
pressuretime = 500                          # Time step duration to hold under pressure before releasing.
max6 = 0                                    # Maximum number of particles with 6 connections, more is better.
mass = 1                                    # Mass, as in Newton's equations of motion.
dt = .005                                   # Timestep size, how (slower but) precise it moves.
ljlimit = 10000                             # Limit to amount of force to apply to motion.
epsilon = 1                                 # Lennard Jones variable, remains 1.
sigma = 1                                   # Lennard Jones variable, remains 1.
timestepnumber = 0                          # For displaying the time step number.
lastsave = 250                              # Prevent too many files when saving best (5 or more per point) symmetry.
filenumber = 0                              # Number to apply to text file pathnames saved to disk for each write.

# Initialize global lists
pos = [[0, 0, 0] for _ in range(pointcount)]        # Positions of points.
vel = [[0, 0, 0] for _ in range(pointcount)]        # Velocities.
acc = [[0, 0, 0] for _ in range(pointcount)]        # Accelerations.
force = [[0, 0, 0] for _ in range(pointcount)]      # Amount of Lennard-Jones force to apply to velocity.
ljsum = [0 for _ in range(pointcount)]              # Clear sum used for color and brightness.
xyc = [[] for _ in range(pointcount)]               # X,Y and color to draw on screen.
connect = [[] for _ in range(pointcount)]           # List of what connects to what, neighbors.

def writefiles():
global filenumber
filenumber += 1
filename1 = 'xyz' + str(filenumber) + 'Points.txt'
filename2 = 'xyz' + str(filenumber) + 'Connect.txt'
print('Saving ', filename1, ' and ', filename2, ' file')  # Then save file to disk, same folder program's in.
csvwrite(pos, filename1)
csvwrite(connect, filename2)
global screenfillcolor
screenfillcolor = White                 # Make screen background flash white to indicate file save.
return dialog

def csvwrite(sourcelist, filename):         # Write Comma Separated Value list of lists [,,] or tuples (,,) to disk.
# For control over file behavior csv or other import module is not used for reading or writing.
with open(filename, 'w') as f:          # Open file to write over, using 'f' as its name:
for line in sourcelist:             # For each line of values in source list:
s = str(line)                   # Convert line to a string containing commas.
f.write(s[1:-1] + "\n")         # Write all but [] or () chars, line feed at end.

# File to read must have no blank lines or comments added to it, comma separated values only.
destination = []                        # Clear destination list lists are appended to.
with open(filename, 'r') as f:          # Open file in read mode using 'f' as its name:
for line in f:                      # For each comma separated line of code in 'f':
s = line.split(',')             # Split line at commas into list named 's'
d = []                          # Empty list named 'd' to fill with elements.
for e in s:                     # For each string element listed in 's' list:
d.append(float(e))          # Append floating point element value to list.
destination.append(d)           # Append list of elements to destination list.
return destination                      # Returns list of lists only = [[,,],[,,],[,,]]

#  Update positions, velocities, accelerations.
def update(points, pos, vel, f, acc, mass, dt):
rmass = 1.0 / mass
if dialog != 'Holding Own Shape':
tosphere()
# Update positions, for 3 dimensions.
for i in range(pointcount):
for j in range(3):
pos[i][j] = pos[i][j] + vel[i][j] * dt + 0.5 * acc[i][j] * dt * dt
if dialog != 'Holding Own Shape':
tosphere()
#  Update velocities.
for i in range(pointcount):
for j in range(3):
vel[i][j] = vel[i][j] + 0.5 * dt * (f[i][j] * rmass + acc[i][j])
vel[i][j] = vel[i][j]*.2
#  Update accelerations.
for i in range(pointcount):
for j in range(3):
acc[i][j] = f[i][j] * rmass
return pos, vel, acc

def forces(points, pos, vel, mass):
#  Compute Lennard-Jones force and energies.
ljt = 0  # Initialize variable for totaling Lennard-Jones force.
ljmax = -100000000                              # Initialize maximum, negative when points attract.
ljmin = 1000000000000000000000000000000000      # Minimum reading goes negative when points attract.
ljsum = [0 for _ in range(points)]              # Clear sum used for color and brightness drawn to screen.
force = [[0,0,0] for _ in range(points)]
rij = [0,0,0]
global connect
connect = [[] for _ in range(points)]  # Positions of points.
distribution = [0 for _ in range(10)]           # list of how many in total have 0,1,2,3,4,5,6,7,8 and 9 connections.
potential = 0.0
for i in range(len(near)):
for fj in range(len(near[i])):
j = near[i][fj]
#  Displacement vector RIJ.
for k in range(3):
rij[k] = pos[i][k] - pos[j][k]
#  Compute D and D2, a distance and a truncated distance.
d = 0.0
for k in range(3):
d = d + rij[k] ** 2
d = np.sqrt(d)
if d < maxconnectspace:
connect[i].append(j)
#  Lennard Jones force.
ljr = d / spacing  # Lennard-Jones radius. At ljr=1 force is 0.
if ljr <= 0:  # Prevent divide by zero error when LJ radius is zero.
lj = 0
else:
lj = 4 * epsilon * ((sigma / ljr) ** 12 - (sigma / ljr) ** 6)
if lj>ljlimit:
lj=ljlimit
ljsum[i] += lj                  # Used for brightness to draw on screen.
#  Add particle J's contribution to the force on particle I.
f = lj
for k in range(3):
force[i][k] = force[i][k] + (rij[k] * f)
# Add force to sum total, for average of all points.
ljt += ljsum[i]
# Save minimum and maximum readings.
if ljsum[i] > ljmax:
ljmax = ljsum[i]
if ljsum[i] < ljmin:
ljmin = ljsum[i]
# Divide down sum total for average.
if nearcount > 0:
ljavr = ljt / nearcount
else:
ljavr = 0
# Save distribution for 1 ro 10 connections, most should have 6.
for i in range(points):
connectioncount = len(connect[i])                   # List length at i same as number of connections.
if connectioncount < len(distribution):             # Normally only 0 to 9 will ever be nonzero.
distribution[connectioncount] += 1              # Increment place in the list, for given connections.
return force, potential, ljavr, ljmax, ljmin, ljsum, distribution

# Cartesian to Spherical coordinates.
def cart2sphere(xyz):
hxy = np.hypot(xyz[0], xyz[1])
el = np.arctan2(xyz[2], hxy)
az = np.arctan2(xyz[1], xyz[0])
return az, el, np.hypot(hxy, xyz[2])

# Spherical to Cartesian coordinates.
def sphere2cart(az, el, ra):
theta = ra * np.cos(el)
x = theta * np.cos(az)
y = theta * np.sin(az)
z = ra * np.sin(el)
return [x, y, z]

def drawaxisline(x, y, z, c):
az, el, r = cart2sphere([x, y, z])                  # Cartesian to Spherical (az, el, r)
x1, y1, z1 = sphere2cart(az + rotation, el, r)      # Spherical to Cartesian.
pygame.draw.line(screen, c, (xc, yc),  (x1+xc, z1+yc), 1)

def update_nearbylist():
near = [[] for _ in range(pointcount)]  # Initialize list of points within cutoff range.
nearcount = 0  # And a counter for keeping track of how many.

for i in range(pointcount):                         # For all the i points
for j in range(pointcount):                     # and for all the j points
if i != j:                                  # Except itself
dx = pos[j][0] - pos[i][0]              # Calculate x,y,z distances.
dy = pos[j][1] - pos[i][1]
dz = pos[j][2] - pos[i][2]
# Use the three distances for calculating an always positive 3D distance.
# This is the very time consuming command, to less occasionally use.
dist = np.sqrt((dx * dx) + (dy * dy) + (dz * dz))
# If distance is within cutoff range then save 'j' point, at this 'i'.
if dist < cutoffdist:
near[i].append(j)
nearcount += 1
print('Computed Nearby List of Points within CutOff Range, ', nearcount, 'distances')
return near, nearcount

def render():
screen.fill(screenfillcolor)
# Precalculate all the rotated X,Y screen locations and their colors.
for i in range(pointcount):
red255 = 0
grn255 = 0
blu255 = 0
az, el, r = cart2sphere(pos[i][0:])    # Cartesian to Spherical (az, el, r)
x, y, z = sphere2cart(az + rotation, el, r)     # Spherical to Cartesian.
# Vary brightness by amount of force, positive increasingly orange, negative blue.
if ljsum[i] > 0:
a1 = (ljsum[i]/ljlimit)* HalfPi
if a1 > HalfPi:
a1 = HalfPi
red255 = np.sin(a1)*255
grn255 = red255*.6
blu255 = 0
else:
red255 = 0
grn255 = red255
blu255 = int(-ljsum[i] * 20)
if blu255 > 255:
blu255 = 255
c = (red255, grn255, blu255)
xyc[i] = [x, y, z, c]               # Save xyz and color.
# To properly draw on top of each other sort the list according to what is furthest away, [1] = by Y axis.
xyc.sort(key=lambda x: x[1])
# Draw all points to screen.
for i in range(pointcount):
# Calculate x,y and radius then draw point to screen buffer.
x, y, z, c = xyc[i]
xs = int(x + xc)
ys = int(z + yc)
if xs > 10000:
xs = 10000
if ys > 10000:
ys = 10000
cr = int(spacing/2)
if y > -spacing:
pygame.draw.circle(screen, c, [xs, ys], cr, 0)
pygame.draw.circle(screen, (255, 255, 255), [xs, ys], cr, 1)
# Draw axis lines.
ra = sphereradius * 1.1                 # Distance of axis line, from 0.
drawaxisline(0, 0, ra, Yellow)
if rotation > ThreeQtrRadian:           # Green x behind blue y?
drawaxisline(ra, 0, 0, Green)
drawaxisline(0, ra, 0, Blue)
if rotation < ThreeQtrRadian:           # Green x on top of blue y?
drawaxisline(ra, 0, 0, Green)
pygame.display.update()                     # Display everything drawn into screen buffer.
fps.tick(60)

def tosphere():
# Keep points on sphere by converting to Spherical then (using sphere radius) to Cartesian again.
for i in range(pointcount):                     # 1D Radial distance between x points.
az, el, r = cart2sphere(pos[i])             # Cartesian to Spherical.
xyz = sphere2cart(az, el, sphereradius)     # Spherical to Cartesian again.
pos[i] = xyz

# Start program.
fnd = 1                                             # Calculate again right away after starting.
dialog = "Squeezing"                                # Set to begin squeeze cycle by changing dialog.
holdcount = pressuretime                            # Time to hold squeeze for.

# Randomly place points anywhere in array list. Points will later become attached to surface of sphere.
for i in range(pointcount):
az = random.uniform(-Pi, Pi)                    # Guess fractional number between -3.14 to +3.14 for azimuth.
el = random.uniform(-Pi, Pi)                    # Guess fractional number between -3.14 to +3.14 for elevation.
xyz = sphere2cart(az, el, sphereradius)         # Spherical to Cartesian.
pos[i] = xyz

# PyGame Initialization
pygame.init()
screen = pygame.display.set_mode(screensize)
xc = int(screensize[0]/2)
yc = int(screensize[1]/2)
pygame.display.set_caption("Vesicle  s key to save,  Pos. Axis: x=green, y=blue, z=yellow")
fps = pygame.time.Clock()
paused = False

while True:
screenfillcolor = Black  # Normally black screen background.
for event in pygame.event.get():
if event.type == pygame.QUIT:
pygame.quit()
sys.exit()
if event.type == pygame.KEYUP:
if event.key == pygame.K_SPACE:
paused = not paused
if event.key == pygame.K_s:
writefiles()
if not paused:
timestepnumber += 1
# Update positions of points, and rotation to drawn on screen.
fnd += 1            # Advance count for next update to nearby within cutoff range list.
if fnd > (cutoffcalctime * (1/(2000/timestepnumber))):              # Updates less often over time.
fnd = 0                                                         # Start new count.
near, nearcount = update_nearbylist()                           # Update nearby, within cutoff range.
pos, vel, acc = update(pointcount, pos, vel, force, acc, mass, dt)  # Update xyz location of each location.
force, potential, ljavr, ljmax, ljmin, ljsum, dn = forces(pointcount, pos, vel, mass)   # Update forces.

# Squeeze by reducing radius of sphere, by average L-J force amount.
moreless = ""                               # Clear optional "More" or "Less" dialog.
if dialog == "Squeezing":
holdcount -= 1                          # Decrement squeezing time counter.
if holdcount < 0:                       # Squeeze time ended?
holdcount = pressuretime            # Start countdown to next expansion cycle.
dialog = "Begin Expansion LJ=0"     # Changing 'dialog' ends Squeezing.
if ljavr < pressuremax and ljmax < ljlimit:
moreless = "More, "                 # Add to dialog as "Squeezing More".
else:
moreless = "Less, "                 # Add to dialog as "Squeezing Less".
else:
dialog = "Squeezing"                # Changing dialog to "Squeezing" ends expansion.
writefiles()                        # Save file, of expanded state.
else:
if ljavr > 0:
else:
print(timestepnumber, "   ", " Connection Distribution", dn, dialog, moreless, "Radius =", int(sphereradius), "  LJaverage =", ljavr, "  LJmax =", ljmax, "  LJmin =", ljmin)

#       fiveOrMore = dn[0] + dn[1] + dn[2] + dn[3] + dn[4]
if dn[3] == 0 and dn[4] == 0 and max6 <= dn[6] and dn[7] == 0 and lastsave + 20 < timestepnumber:
max6 = dn[6]
#            if fiveOrMore == 0 and ljavr < 5 and lastsave+20 < timestepnumber:
#                if fiveOrMore == 0 and ljavr < 5 and lastsave + 20 < timestepnumber:
lastsave = timestepnumber
writefiles()

# Print everything to screen.
render()
rotation += rotationstep                    # Add to screen rotation amount.
if rotation > Radian:                       # Keep in range from 0 to 1 Radian
rotation = 0``````

Zip file containing VesicleMaker.py and VesicleFileViewer.py is at:

Vesicle maker screen flashes when saving files with good 5 and 6 symmetry and periodically every 500 timesteps or so, to meanwhile scan through full series using left and right arrow keys. Two example files showing start and after 676 timesteps will be written over after starting the Maker program.

Nice looking articulation files are now being generated. Imperfections that may still exist (possibly due to needing to add or subtract particles to make perefectly symmetrical all around) will be part of the test, question is how well the wave signal travels through less than perfect symmetry.

Since this applies to living cell soma and hillock membrane signal propagation (and HTM Theory considers each “cell” in the memory to be a separate predictive unit) it may be the basis for a novel molecular biology level 360 degree spatial “HTM neuron” that is controlled by usual passing, reflecting or (do nothing) blocking propagation, along the path towards axon activation. Already have a vesicle for soma, and can add bulbous hillock then string of particles for axon outputs and dendrites. Molecular biologists might be able to from there help get something working.

As in drawing a map of physical space where there are things to bump into and avoid complex spatial predictions were intuitive. Path it needs to take at the time becomes pointed out as direction to go when there.

For me wave action made complex behaviors that could never be “coded in” some other way come to life, and may again, by including cell membrane signal propagation even though it sounds harder that way. What you know from modeling neurons as you do this way become clues as to how cells radially control signal flow along their length.

It’s again likely not exactly HTM looking, but this is at least something new that works with HTM and Thousand Brains Theory, to help us think forward on into biology. I hope that with my having to get accustomed to Python ways and needing lots of long experiments this turned into way more of a project than I origionally planned this would make a last minute gift for under your Christmas tree. I now only need to light one up, with some wave action!

2 Likes

Signal propagation works better than I expected! Code is at: https://sites.google.com/site/intelligencedesignlab/home/BrainBall.zip

Even with way less than perfect vesicle symmetry the waves are soon cancelling out real nice on the opposite side of the sphere.

There is no longer a file reader program. The maker program shows front and back sides, both ways, with signals on top.

If there is a “xyzStartup.txt” file (saved by pressing “s” key) then it uses that, otherwise program starts from scratch on a new one.

When an area starts bouncing around enough to break and reform connections the result is overactivity resembling an epileptic seizure. Since it is of interest to see what happens when connections go to the wrong places I left that in, instead of preventing by silencing cells that are breaking then reforming connections. Program will clear all outputs after becoming too active, but shows enough activity to get an idea what the possible self-oscillation patterns look like. We then get to see complex wave patterns that can on purpose also be generated, for some purpose, like this spiral:

The input pattern each receives normally indicates where on the sphere the wave originated from. On the direct opposite side of the sphere all inputs become 1 then cancel after negating input, to derive output. Input bits 111111 become output bits 000000, no output signal. It’s one of the things (assuming they are intelligent enough) early cells could possibly do for fun “because they can” then later use for polled navigational planning or other purposes.

Full code is:

``````# Propagate signals around a developing "brain ball" vesicle.  Gary Gaulin, 1/1/2020
# Lennard Jones Potential approximates stem cell migration on a sphere to equally space themselves apart.
# If there is a "xyzStartup.txt" file (saved by pressing "s" key) then program starts with that, otherwise from scratch.
# There is also a "Nice-xyzStartup.txt" from earlier benchmark to use to improve upon, by tweaking physics related code.
# It is not necessary to achieve perfect symmetry, waves will still travel fine, but symmetry can be improved.
# This implementation is for also showing what happens when connections are breaking and signals go to wrong places,
# which can in the timestepsignals function be prevented, when connections change from one timestep to next.
# Keyboard space bar pauses program, then resumes.
# Thin 3D Axis lines are drawn in upper left view showing: x=green, y=blue, z=yellow
# https://discourse.numenta.org/t/python-program-using-lennard-jones-force-to-approximate-cortical-stem-cell-close-packing-symmetry/6456

import random
import sys
import numpy as np
import pygame

Pi = 3.14159265358979323846264338327950288412
HalfPi = Pi / 2
Black = (0, 0, 0)
Brown = (55, 0, 0)
Red = (255, 0, 0)
Orange = (155, 100, 50)
Yellow = (220, 220, 0)
Green = (0, 200, 0)
Blue = (0, 0, 150)
Violet = (255, 0, 255)
Gray = (120, 120, 120)
White = (255, 255, 255)
# From electronics: 0=Black, 1=Brown, 2=Red, 3=Orange, 4=Yellow, 5=Green, 6=Blue, 7=Violet, 8=Gray, 9=White
ElectronicColorCode = Black, Brown, Red, Orange, Yellow, Green, Blue, Violet, Gray, White
# Colors used to draw connections that are either 0, or 1 for action potential.
SignalColor = Black, White

screensquare = 360                          # Size for each of 4 squares.
screensize = Width, Height = (screensquare * 2, screensquare * 2)    # PyGame window size, to show 2*2=4 views.
xc = int(screensize[0]/4)                   # Center x of 1 screen.
yc = int(screensize[1]/4)                   # Center y of 1 screen.
timestepnumber = 0                          # For displaying the time step number.

def csvwrite(sourcelist, filename):         # Write Comma Separated Value list of lists [,,] or tuples (,,) to disk.
# For control over file behavior csv or other import module is not used for reading or writing.
with open(filename, 'w') as f:          # Open file to write over, using 'f' as its name:
for line in sourcelist:             # For each line of values in source list:
s = str(line)                   # Convert line to a string containing commas.
f.write(s[1:-1] + "\n")         # Write all but [] or () chars, line feed at end.

# File to read must have no blank lines or comments added to it, comma separated values only.
destination = []                        # Clear destination list lists are appended to.
with open(filename, 'r') as f:          # Open file in read mode using 'f' as its name:
for line in f:                      # For each comma separated line of code in 'f':
s = line.rstrip('\n')
d = []                          # Empty list named 'd' to fill with elements.
if len(s) > 0:
s = s.split(',')            # Split line at commas into list named 's'
for e in s:                 # For each string element listed in 's' list:
d.append(float(e))      # Append floating point element value to list.
destination.append(d)           # Append list of elements to destination list.
return destination                      # Returns list of lists only = [[,,],[,,],[,,]]

#  Update positions, velocities, accelerations.
def update(points, pos, vel, f, acc, mass, dt):
rmass = 1.0 / mass
if dialog != 'Holding Own Shape':
tosphere()
# Update positions, for 3 dimensions.
for i in range(pointcount):
for j in range(3):
pos[i][j] = pos[i][j] + vel[i][j] * dt + 0.5 * acc[i][j] * dt * dt
if dialog != 'Holding Own Shape':
tosphere()
#  Update velocities.
for i in range(pointcount):
for j in range(3):
vel[i][j] = vel[i][j] + 0.5 * dt * (f[i][j] * rmass + acc[i][j])
vel[i][j] = vel[i][j]*.2
#  Update accelerations.
for i in range(pointcount):
for j in range(3):
acc[i][j] = f[i][j] * rmass
return pos, vel, acc

def forces(points, pos, vel, mass):
#  Compute Lennard-Jones force and energies.
ljt = 0  # Initialize variable for totaling Lennard-Jones force.
ljmax = -100000000                              # Initialize maximum, negative when points attract.
ljmin = 1000000000                              # Minimum reading goes negative when points attract.
ljsum = [0 for _ in range(points)]              # Clear sum used for color and brightness drawn to screen.
force = [[0, 0, 0] for _ in range(points)]
rij = [0, 0, 0]
global connect
connect = [[] for _ in range(points)]           # Articulation list for what connects to where.
distribution = [0 for _ in range(10)]           # list of how many total have 0,1,2,3,4,5,6,7,8 and 9 connections.
lj = 0
for i in range(len(near)):
for fj in range(len(near[i])):
j = near[i][fj]
#  Displacement vector RIJ.
for k in range(3):
rij[k] = pos[i][k] - pos[j][k]
#  Compute D and D2, a distance and a truncated distance.
d = 0.0
for k in range(3):
d = d + rij[k] ** 2
d = np.sqrt(d)
connect[i].append(j)
#  Lennard Jones force.
ljr = d / spacing  # Lennard-Jones radius. At ljr=1 force is 0.
if ljr <= 0:  # Prevent divide by zero error when LJ radius is zero.
lj = 0
else:
lj = 4 * epsilon * ((sigma / ljr) ** 12 - (sigma / ljr) ** 6)
if lj > ljlimit:
lj = ljlimit
lj = lj / ljdivide
ljsum[i] += lj                  # Used for brightness to draw on screen.
#  Add particle J's contribution to the force on particle I.
for k in range(3):
force[i][k] = force[i][k] + (rij[k] * lj)
# Add force to sum total, for average of all points.
ljt += ljsum[i]
# Save minimum and maximum readings.
if lj > ljmax:
ljmax = lj
if lj < ljmin:
ljmin = lj
# Divide down sum total for average.
if nearcount > 0:
ljavr = ljt / nearcount
else:
ljavr = 0
# Save distribution for 0 to 9 connections, most should evetually have 6.
for i in range(points):
connectioncount = len(connect[i])                   # List length at i same as number of connections.
if connectioncount < len(distribution):             # Normally only 0 to 9 will ever be nonzero.
distribution[connectioncount] += 1              # Increment place in the list, for given connections.
return force, ljavr, ljmax, ljmin, ljsum, distribution

# Cartesian to Spherical coordinates.
def cart2sphere(xyz):
hxy = np.hypot(xyz[0], xyz[1])
el = np.arctan2(xyz[2], hxy)
az = np.arctan2(xyz[1], xyz[0])
return az, el, np.hypot(hxy, xyz[2])

# Spherical to Cartesian coordinates.
def sphere2cart(az, el, ra):
theta = ra * np.cos(el)
x = theta * np.cos(az)
y = theta * np.sin(az)
z = ra * np.sin(el)
return [x, y, z]

def drawaxisline(x, y, z, c):
az, el, r = cart2sphere([x, y, z])                  # Cartesian to Spherical (az, el, r)
x1, y1, z1 = sphere2cart(az + rotation, el, r)      # Spherical to Cartesian.
pygame.draw.line(screen, c, (xc, yc),  (x1+xc, z1+yc), 1)

def update_nearbylist():
near = [[] for _ in range(pointcount)]  # Initialize list of points within cutoff range.
nearcount = 0  # And a counter for keeping track of how many.
for i in range(pointcount):                         # For all the i points
for j in range(pointcount):                     # and for all the j points
if i != j:                                  # Except itself
dx = pos[j][0] - pos[i][0]              # Calculate x,y,z distances.
dy = pos[j][1] - pos[i][1]
dz = pos[j][2] - pos[i][2]
# Use the three distances for calculating an always positive 3D distance.
# This is the very time consuming command, to less occasionally use.
dist = np.sqrt((dx * dx) + (dy * dy) + (dz * dz))
# If distance is within cutoff range then save 'j' point, at this 'i'.
if dist < cutoffdist:
near[i].append(j)
nearcount += 1
print('Computed Nearby List of Points within CutOff Range, ', nearcount, 'distances')
return near, nearcount

def render():
# Precalculate all the rotated X,Y screen locations and their colors.
for i in range(pointcount):
az, el, r = cart2sphere(pos[i][0:])             # Cartesian to Spherical (az, el, r)
x, y, z = sphere2cart(az + rotation, el, r)     # Spherical to Cartesian.
# Vary brightness by amount of force, positive increasingly orange, negative blue.
if ljsum[i] > 0:
a1 = (ljsum[i]/ljlimit) * ljdivide * HalfPi
if a1 > HalfPi:
a1 = HalfPi
red255 = np.sin(a1)*255
grn255 = red255*.6
blu255 = red255*.2
else:
red255 = 0
grn255 = red255
blu255 = int((-ljsum[i]/10) * ljdivide * 255)
if blu255 > 255:
blu255 = 255
c1 = (red255, grn255, blu255)
n = len(connect[i])
if n <= 9:
c2 = ElectronicColorCode[n]
else:
c2 = White
xyzc[i] = [x, y, z, c1, c2, i]                  # Save xyz and color.
unsortedxyzc[i] = xyzc[i]
# To properly draw on top of each other sort the list according to what is furthest away, [1] = by Y axis.
xyzc.sort(key=lambda x: x[1])
cr = int(spacing / 2)
# Draw two screens on left shoiwing front view.
for i in range(pointcount):
x, y, z, c1, c2, n = xyzc[i]                    # Normal i direction, nearest is drawn on top of far away.
xs = int(x + xc)
ys = int(z + yc)
if y > -spacing:                        # For clarity only nearest/visible side of sphere is shown.
pygame.draw.circle(screen, c1, [xs, ys], cr, 0)
pygame.draw.circle(screen, (255, 255, 255), [xs, ys], cr, 1)
pygame.draw.circle(screen, c2, [xs, ys + screensquare], cr, 0)
pygame.draw.circle(screen, (255, 255, 255), [xs, ys + screensquare], cr, 1)
# Draw two screens to the right showing back views.
for i in range(pointcount):
x, y, z, c1, c2, n = xyzc[pointcount-1-i]       # Reversed i direction, far away is instead drawn on top.
xs = int(-x + xc)                               # Viewed from back: what was on left is to right, negative x.
ys = int(z + yc)                                # Z axis remains same up/down orientation.
if y < spacing:                                 # For clarity only nearest/visible side of sphere is shown.
pygame.draw.circle(screen, c1, [xs + screensquare, ys], cr, 0)
pygame.draw.circle(screen, (255, 255, 255), [xs + screensquare, ys], cr, 1)
pygame.draw.circle(screen, c2, [xs + screensquare, ys + screensquare], cr, 0)
pygame.draw.circle(screen, (255, 255, 255), [xs + screensquare, ys + screensquare], cr, 1)
# Draw signal propagation lines on top of 2 front views.
for i in range(pointcount):
x, y, z, c1, c2, n = xyzc[i]
if y > -spacing:                         # For clarity only nearest/visible side of sphere is shown.
xs1 = int(x + xc)
ys1 = int(z + yc)
bits = outs[n]
for j in range(len(connect[n])):
x2, y2, z2, c1, c2, n2 = unsortedxyzc[connect[n][j]]
drawsignalline(xs1, ys1, x2 + xc, z2 + yc, bits[j], 0, screensquare)
drawsignalline(xs1, ys1, x2 + xc, z2 + yc, bits[j], 0, 0)
# Draw signal propagation lines on top of 2 back views.
for i in range(pointcount):
x, y, z, c1, c2, n = xyzc[pointcount - 1 - i]
if y < spacing:                          # For clarity only nearest/visible side of sphere is shown.
xs1 = -x + xc
ys1 = z + yc
bits = outs[n]
for j in range(len(connect[n])):
x2, y2, z2, c1, c2, n2 = unsortedxyzc[connect[n][j]]
drawsignalline(xs1, ys1, -x2 + xc, z2 + yc, bits[j], screensquare, screensquare)
drawsignalline(xs1, ys1, -x2 + xc, z2 + yc, bits[j], screensquare, 0)
# Draw thin axis lines in upper left view.
ra = sphereradius * 1.1                     # Distance of axis line, from 0.
drawaxisline(0, 0, ra, Yellow)
if rotation > ThreeQtrRadian:               # Green x behind blue y?
drawaxisline(ra, 0, 0, Green)
drawaxisline(0, ra, 0, Blue)
if rotation < ThreeQtrRadian:               # Green x on top of blue y?
drawaxisline(ra, 0, 0, Green)
ra = sphereradius * 1.1                     # Distance of axis line, from 0.
drawaxisline(0, 0, ra, Yellow)
if rotation > ThreeQtrRadian:               # Green x behind blue y?
drawaxisline(ra, 0, 0, Green)
# Draw lines to divide screen into four.
pygame.draw.line(screen, Gray, (screensquare - 1, 0), (screensquare - 1, Width), 1)
pygame.draw.line(screen, White, (0, screensquare - 1), (Height, screensquare - 1), 1)
# Display everything to screen.
pygame.display.update()                     # Display everything drawn, into screen buffer.
fps.tick(60)                                # Speed is at most 60 frames per second.

def drawsignalline(x1, y1, x2, y2, b, xleft, yleft):
thick = int((b * spacing / 4) + 1)
pygame.draw.line(screen, SignalColor[b], [xleft+x1, yleft+y1], [xleft+x2, yleft+y2], thick)

# Keep points on sphere by converting to Spherical then (using sphere radius) to Cartesian again.
def tosphere():
for i in range(pointcount):                     # 1D Radial distance between x points.
az, el, r = cart2sphere(pos[i])             # Cartesian to Spherical.
xyz = sphere2cart(az, el, sphereradius)     # Spherical to Cartesian again.
pos[i] = xyz

def timestepsignals():
signalactivity = 0
global outs
ins = [[] for _ in range(pointcount)]           # Clear Input list output bits from neighbors are gathered into.
for i in range(pointcount):                     # For each 'i' point in vesicle, each with connections to 'j's
for j in range(len(connect[i])):            # For each of the (usually 6) 'j' connections to other cells,
if j < len(outs[i]):                    # If output connection still exists in connection list then
b = outs[i][j]                      # bit to save is from output list.
else:                                   # Else
b = 0                               # to prevent program error append with 0 so any bit fills void.
ins[connect[i][j]].append(b)            # Append bit from 'j' output to corresponding input line at 'i'
# Use gathered Input bits to set corresponding Output list bits.
outs = [[] for _ in range(pointcount)]          # Clear Output list that new Output bits will be placed into.
for i in range(pointcount):                     # For each 'i' point in vesicle, each with connections to 'j's
signalactivity += sum(ins[i])               # Signal activity equals total number of bits set.
# If any of the Input bits are a 1 then keep wave going by Outputting the opposite bit pattern.
if sum(ins[i]) > 0:                         # If there are any bits set then
outs[i] = [int(not b) for b in ins[i]]  # all 1's in output list become 0, and all 0 becomes 1, negates.
else:                                       # Else
outs[i] = ins[i]                        # Output bits remain all 0.
# If no signal activity then start wave by setting all output connections of a randomly selected place to 1.
if signalactivity == 0:                         # If no signal acctivity at all then
i = random.randint(0, pointcount-1)         # pick one of the cells at random.
outs[i] = [1] * len(connect[i])             # Set all its outputs to 1, to fire action potentials.
signalactivity = len(connect[i])            # Signal activity equals number of bits that were just set.
# If too much uncontrolled signal activity then clear all outputs.
if signalactivity > pointcount*2:               # If more than average two connections per cell then
for i in range(pointcount):                 # for all cells
outs[i] = [0] * len(outs[i])            # fill list with 0's, same length as there are out connections.
print("All outs set 0 to end Epileptic type overactivity from broken/rewiring connections.")

#######################################################################################################################
# Start program.
try:
pointcount = len(pos)
except IOError:
# Randomly place points anywhere in array list. Points will later become attached to surface of sphere.
# Example of 392 is from dimple count of a 1950's Ben Hogan-5 golfball.
# It has unusually uniform symmetry, as though well planned to close pack dimples neatly around a sphere.
pointcount = 392
pos = [[0, 0, 0] for _ in range(pointcount)]    # Positions of points.
for i in range(pointcount):
az = random.uniform(-Pi, Pi)                # Guess fractional number between -3.14 to +3.14 for azimuth.
el = random.uniform(-1.2, 1.2)              # Guess fractional number between -3.14 to +3.14 for elevation.
xyz = sphere2cart(az, el, sphereradius)     # Spherical to Cartesian.
pos[i] = xyz                                # Add x,y,z to list of positions.

sar = sas/pointcount                        # Divide for area of one Lennard-Jones radius, between points.
rbp = np.sqrt(sar/Pi)*2                     # Radius between points, for given area.
spacing = rbp-((1.2/pointcount)*rbp)        # Radius between points, spacing to also fit screen.
connectdist = spacing*1.2                   # Distance two points are touching, bumping each other.
cutoffdist = spacing*2.5                    # Distance two points are too far away to influence each other.
cutoffcalctime = 40                         # Time steps between updates of i,j points within cut-off range.
rotationstep = .02                          # Amount of rotation per timestep.
rotation = 0                                # Initialize rotation angle, for drawing all points on screen.
max6 = 0                                    # Maximum numbrer of particles with 6 connections, more is better.
mass = 1                                    # Mass, as in Newton's equations of motion.
dt = .3                                     # Timestep size, how (slower but) precise it moves.
epsilon = 1                                 # Lennard Jones variable, remains 1.
sigma = 1                                   # Lennard Jones variable, remains 1.
pressurestep = spacing/250                  # Amount of radius change per time step, in screen pixels.
ljlimit = 700                               # Limit to amount of force to apply to motion.
ljdivide = 300                              # Limit to amount of force to apply to motion.
ljmaxforce = ljlimit/ljdivide / 7
ljavrmax = .023                             # Maximum amount of average Lennard-Jones force to apply.
signalactivity = 0
lastsave = 500                              # Prevent too many files when saving best (5 or more per point) symmetry.
fnd = 1                                     # Calculate again right away after starting.
dialog = "Squeezing"                        # Set to begin squeeze cycle by changing dialog.
moreless = ""                               # Clear "More" or "Less" dialog.

# Initialize global lists
vel = [[0, 0, 0] for _ in range(pointcount)]        # Velocities.
acc = [[0, 0, 0] for _ in range(pointcount)]        # Accelerations.
force = [[0, 0, 0] for _ in range(pointcount)]      # Amount of Lennard-Jones force to apply to velocity.
ljsum = [0 for _ in range(pointcount)]              # Clear sum used for color and brightness.
xyzc = [[] for _ in range(pointcount)]              # X,Y and color to draw on screen.
unsortedxyzc = [[] for _ in range(pointcount)]      # Saved before sorting.
connect = [[0] for _ in range(pointcount)]          # List with list of neighboring point numbers to connect to.
connectto = []                                      # List with list of neighboring point numbers to connect to.
outs = [[] for _ in range(pointcount)]              # Output (action potentials) bits, to each 1 bit connection.
ins = [[] for _ in range(pointcount)]               # Input (action potential) bits, from each 1 bit connection.

# PyGame Initialization
pygame.init()
screen = pygame.display.set_mode(screensize)
pygame.display.set_caption("Vesicle Maker Program. Gary Gaulin, 2020")
fps = pygame.time.Clock()
paused = False

while True:
for event in pygame.event.get():
if event.type == pygame.QUIT:
pygame.quit()
sys.exit()
if event.type == pygame.KEYUP:
if event.key == pygame.K_SPACE:
paused = not paused
if event.key == pygame.K_s:
csvwrite(pos, 'xyzStartup.txt')
screen.fill(White)  # Make screen background flash white to indicate file save.
if not paused:
timestepnumber += 1
# Update positions of points, and rotation to drawn on screen.
fnd += 1            # Advance count for next update to nearby within cutoff range list.
if fnd > (cutoffcalctime * (1/(2000/timestepnumber))) or fnd > 100:     # Updates less often over time.
fnd = 0                                                             # Start new count.
near, nearcount = update_nearbylist()                               # Update nearby, within cutoff range.
pos, vel, acc = update(pointcount, pos, vel, force, acc, mass, dt)      # Update xyz location of each location.
force, ljavr, ljmax, ljmin, ljsum, dn = forces(pointcount, pos, vel, mass)  # Forces and connections
timestepsignals()
# Squeeze by reducing radius of sphere, by average L-J force amount.
if dialog == "Squeezing":
if ljavr < ljavrmax and ljmax < ljmaxforce:
moreless = "More, "  # Add to dialog as "Squeezing More".
else:
moreless = "Less, "  # Add to dialog as "Squeezing Less".
print(timestepnumber, "  distribution 0-9=", dn, dialog, moreless, "  R=", sphereradius, "  LJavr=", ljavr, "  LJmax =", ljmax, "  LJmin =", ljmin)
# Print everything to screen.
render()                            # Draw to screen using currently set background color.
screen.fill(Black)                  # Normally black background, gets set to render in white by file save.
rotation += rotationstep            # Add to screen rotation amount.
if rotation > Radian:               # Keep in range from 0 to 1 Radian
rotation = 0``````

The timestep code for starting waves and signal propagation was reduced down to only this:

``````def timestepsignals():
signalactivity = 0
global outs
ins = [[] for _ in range(pointcount)]           # Clear Input list output bits from neighbors are gathered into.
for i in range(pointcount):                     # For each 'i' point in vesicle, each with connections to 'j's
for j in range(len(connect[i])):            # For each of the (usually 6) 'j' connections to other cells,
if j < len(outs[i]):                    # If output connection still exists in connection list then
b = outs[i][j]                      # bit to save is from output list.
else:                                   # Else
b = 0                               # to prevent program error append with 0 so any bit fills void.
ins[connect[i][j]].append(b)            # Append bit from 'j' output to corresponding input line at 'i'
# Use gathered Input bits to set corresponding Output list bits.
outs = [[] for _ in range(pointcount)]          # Clear Output list that new Output bits will be placed into.
for i in range(pointcount):                     # For each 'i' point in vesicle, each with connections to 'j's
signalactivity += sum(ins[i])               # Signal activity equals total number of bits set.
# If any of the Input bits are a 1 then keep wave going by Outputting the opposite bit pattern.
if sum(ins[i]) > 0:                         # If there are any bits set then
outs[i] = [int(not b) for b in ins[i]]  # all 1's in output list become 0, and all 0 becomes 1, negates.
else:                                       # Else
outs[i] = ins[i]                        # Output bits remain all 0.
# If no signal activity then start wave by setting all output connections of a randomly selected place to 1.
if signalactivity == 0:                         # If no signal acctivity at all then
i = random.randint(0, pointcount-1)         # pick one of the cells at random.
outs[i] = [1] * len(connect[i])             # Set all its outputs to 1, to fire action potentials.
signalactivity = len(connect[i])            # Signal activity equals number of bits that were just set.
# If too much uncontrolled signal activity then clear all outputs.
if signalactivity > pointcount*2:               # If more than average two connections per cell then
for i in range(pointcount):                 # for all cells
outs[i] = [0] * len(outs[i])            # fill list with 0's, same length as there are out connections.
print("All outs set 0 to end Epileptic type overactivity from broken/rewiring connections.")``````

Now it’s possible to see the signal dynamics early 1 membrane cell colonies would have to work with, when they only had one input and output connection to share with each of their neighbors. Such a thing would not fossilize like bone but from what I read microscopic traces suggest multicellular animals began as mostly undifferentiated blobs. A common example (among many) the volvox includes traits of plants yet has animal-like internal reproduction:

All this in turn brings us back to the old HTM challenge Matt started that sent me looking for a way to not need hexagonal (or even diagonal) math only good for flat surfaces anyway and has so many methods it could take months or never to settle “on the right one”, which had me back to forming vesicles as a possible way to in the long run save time getting something to on its own explore and recall a simple pygame environment.

For this forum the question is now how to make each cell a predictive “thousand brains” element where each has an “eye spot” for through-straw view of the world around them that moves as they do, so they end up together controlling where they are looking, move in a way all get a good look around them in all or preferred directions. From there scale up to roundish human cortices.

Whatever it takes to get a brain-ball to explore an environment can be a significant clue to the underlying cognitive biology of the most complex brains. Where everything is working right the system goes wherever most of the cells want to go, on its own explores attracting locations with a single cilia hair that beats in response to network wave flow, each cell has some control over by propagating as usual or not. In any case this project led back to the old project Matt started, with at least a novel platform for a 392 sensory unit virtual critter that requires no code to randomly move it around or complex visual and other cortical areas to worry about.

I did not plan when this final stage of the long project would be completed but Python made it easy to work with this kind of bit network, so I ahead of schedule had the vesicle lit up with signal activity, to help celebrate in a New Year with. For at least myself it’s a sigh of relief to have all that behind me, and can now focus on the next step from here.

5 Likes