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
Radian = Pi * 2
ThreeQtrRadian = Radian * 3/4

# 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     #
# points = 32     #
points = 392                        # Total number of (cells, molecules, etc.) points.
screensize = Width, Height = (640, 640)
sr = Height*.38                     # Screen radius, with added space around sphere.
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)
            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():
    spherewidth = sphereradius * 2
# 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., c, [xs, ys], cr, 0), (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.

# 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.
sphereradius = spheremaxradius                  # Start with maximum size sphere.
dialog = "Squeezing"
holdcount = pressuretime

# PyGame Initialization
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:
        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()
        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:
                sphereradius -= pressurestep        # Decrease sphere radius, increase pressure.
                moreless = "More, "                 # Add to dialog as "Squeezing More".
                sphereradius += pressurestep        # Increase sphere radius, decrease pressure.
                moreless = "Less, "                 # Add to dialog as "Squeezing Less".
            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.
                dialog = "Expanding"                # Otherwise use "Expanding" for dialog instead.
                sphereradius += pressurestep        # Increase sphere radius, decreases pressure.
        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.



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.