The post was originally published on CodePen.

This is a pen that's been on my mind for a long time, ever since I started following the work of Anders Hoff that deals with natural growth patterns and emergent structures. I was fascinated by how ordered structures and interesting patterns emerge when simulating simple physical and biological rules. It didn't really seem intuitive how this worked, so I decided to explore the rules that would generate something like the images Anders generated in his own work.

Deconstructing systems

When I started thinking about how to implement such a system, I started off with a list of axioms that would help me model such a system. I came up with the following:

1. An individual cell has no knowledge of the world outside of its immediate neighbourhood.
2. All cells are identical.
3. There is no guiding force governing cell organization other than Newtonian physics and biological growth.

With these axioms in place, I started thinking about how multicellular organisms work in terms of growth. Every organism starts off with a single cell which then divides into two, held together by binding proteins embedded in each of the cells' membrane. For simplicity, let's model cell division by creating a new cell that is attached to its parent by a spring, which represents the connection between the two new cells:

    o  ==(division)==>  o~o   ==(physics)==>  o—o
  (cell)           (the new pair)       (the spring relaxes,
                                        separating the cells)

In this diagram, a squiggly line represents a spring under tension, which pushes the cells apart after relaxation. Sounds simple enough. But what about subsequent divisions? Let's take a look at the next stage of development:

    o—o   ==(division)==>   o~[o~o]
      ^(this cell divides)     ^(the new cell appears)

A subsequent division here causes tension in both springs, which will cause the cells to push against each other to resolve the tension. If we use our intuition to model the physics, we expect the tensions to resolve in one of two ways:

                    o
                   / \
   o~[o~o]   =>   o   o   (when the new cell is a little bit
                                    above the line)

   o~[o~o]   =>  o—o—o    (when the new cell is exactly along
                          the line between the original cells)

NOTE: Here we assume that the spring's direction can freely change, as proteins can flow freely along the membrane in any direction.

So far, it seems like a sensible description of growth for a line of cells. We do have to set a couple of physical rules to the system to model it sensibly:

1. Growth happens in a viscous medium which preventing cells from gliding across space without friction.
2. Cells cannot occupy the same space at the same time, so cells not connected will still repel each other if too close.

Putting it all together

So let's recap the rules:

1. Cells can be connected by springs, tending to keep connected cells at a fixed distance.
2. A cell can divide, which creates a pair of interconnected cells attached to themselves and their original neighbours.
3. Cells are contained in a viscous medium.

We now have all the building blocks we need to build up the models and equations that will govern the evolution of this system.

The cell

Let's start by defining all the properties a cell has:

class Particle {
  constructor(position) {
    this.position = position;
    this.force = new Vec2d(0, 0);
  }
}

The position and force properties are both instances of Vec2d which I use for vector math. We only need those two because we will assume that the medium the cell is embedded in is so viscous that the velocity of the particle drops to 0 as soon as the force ceases to act on the cell (think honey or molasses). During a simulation step, we sum up all the forces acting on the cell, move the particle in the direction of the force and then reset the force property back to 0 for the next step.

The world

We also want to create an environment in which the cells live:

class Simulation {
  constructor(particles) {
    this.particles = particles;
  }
  // ...

Here, we'll keep track of the particles contained in the world. And, to make things simple, we'll assume that particles next to each other in the array are also connected by a spring. That way, when we insert an element into the middle of the array to simulate cell division, the cell will implicitly be connected to its neighbouring cells:

[a,b,c,d] ==> a—b—c—d

But let's make things interesting and assume that cells connect in a circular structure and the last element in the array is attached to the first one:

              a—b
[a,b,c,d] ==> | |
              d—c

So this means that for every element n we have two neighbours n-1 and n+1 connected by springs. If the value exceeds the range [0, particles.length - 1] we wrap the values so they fit the array.

Constants

We define a couple of static constants on the simulation that defines a couple of important values:

  static get NODE_DISTANCE()      { return 2; }
  static get DISTANCE_CUTOFF()    { return this.NODE_DISTANCE * 8; }
  static get SPRING_COEFFICIENT() { return 0.02; }
  static get REPULSION_FORCE()    { return 4; }

Let's go over these constants to get a better understanding about what they signify:

• NODE_DISTANCE: represents the distance (in pixels) between cells that would result in a completely relaxed spring
• SPRING_COEFFICIENT: represents the spring's stiffness; the higher the number, the more force is exerted when the spring is stretched or compressed
• REPULSION_FORCE: the force of repulsion between individual cells, regardless if the're connected or not
• DISTANCE_CUTOFF: represents the maximum distance (in pixels) the cells repel each other; this is to prevent the repulsion from adding up significantly when large numbers of particles are generated.

Of course, these values were chosen completely arbitrarily and through trial and error; you might want to experiement with different values to see how they affect the final shape of the structure.

Calculating forces

We now need to define two types of forces that act on pairs of cells: the spring force and the mutual repulsion.

  static getSpringForce(a, b) {
    let f = Vec2d.sub(a.position, b.position),
        m = f.mag;
    return f.norm().scale((this.NODE_DISTANCE - m) * this.SPRING_COEFFICIENT);
  }

The first function takes two cells and calculates the force that results from the spring acting on the first particle. The equation is simply Hooke's law: F = k(l - m), where l is the relaxed length of the spring (NODE_DISTANCE) and m is the distance between the cells. To calculate the force on the second cell, we just negate the force as prescribed by Newton's Third Law without having to recalculate it.

  static getRepulsionForce(a, b) {
    let f = Vec2d.sub(a.position, b.position),
        m = f.mag;
    return m > this.DISTANCE_CUTOFF ? f.reset(0, 0) : f.norm().scale(this.REPULSION_FORCE / (m * m));
  }

For the repulsion force, we calculate the distance between the cells (m) and use the Inverse-square law to create a force that repels cells from each other: F = k / m^2 (k == REPULSION_FORCE). Because of the inverse square falloff, the cell's influence on others declines quickly with distance, but is very intense when the distance between cells decreases towards 0.

Applying forces and moving particles

Let's move on to going through the particle list and calculating the forces. We make a step(dT) function that will allow us to advance time in the world we've constructed.

Let's start with the beginning:

  step(dT) {
    let t = dT / 4;

Here, we just scale down the time step (passed into the simulation as milliseconds). This can be tuned up or down to slow or speed up the simulation, but if the step is too large, it can lead to instabilities in the simulation. (For a good explanation of errors, see numerical stability in Euler's method.)

    // process mutual repulsion
    for (let i = 0; i < this.particles.length - 1; i++) {
      for (let j = i + 1; j < this.particles.length; j++) {
        // mutual deflection force
        let a = this.particles[i], b = this.particles[j];
        let f = Simulation.getRepulsionForce(a, b);
        a.force.add(f);
        b.force.sub(f);
      }
    }

First, we calculate the repulsion forces between all pairs of cells. The force that getRepulsionForce(a, b) returns is the force acting on particle a, so we add that to the force collector a.force. Due to Newton's Second Law, we must subtract that same force from b.force. This optimization allows us to cut the number of iterations in half, so that this part runs in O(n^2/2) time instead of O(n^2).

    // process linkage
    for (let i = 0; i < this.particles.length; i++) {
      let a = this.particles[i], b = this.particles[(i + 1) % this.particles.length];
      let f = Simulation.getSpringForce(a, b);
      a.force.add(f);
      b.force.sub(f);
    }

This part goes through the chain of elements starting with the first element and calculating the force between it and the next one in the chain. For the last element, it assumes connection with the first one in the array (hence why (i + 1) % this.particles.length). Again, Second Law requires us to add the force to the first cell and subtract from the second.

    // apply movement
    this.particles.forEach(particle => {
      let f = particle.force.clone().scale(t);
      particle.position.add(f);
      particle.force.reset(0, 0);
    });
  }

Finally, we do a bit (okay, A LOT) of cheating: instead of calculating the drag, friction, forces, and integrating twice over the time step to calculate the new position, we instead scale the sum of forces on the particle by the time step and reset the cell's force collector to 0. This simulates the viscosity of the medium while still allowing the cells to move an acceptable distance. (The approach is similar to Euler's method of integration, with the same tradeoff between accuracy and simplicity, and the reason why we want to keep the step size small.)

Initial state

We decided to link the ends of the string of cells into a circular shape, so we need to generate the initial configuration of cells within a circle. We set up an instance of the simulation with a number of particles and a given radius:

// set up the initial world, just start with a circle of 40 particles
var simulation = (() => {
  let n = 40, particles = [], r = 40;
  for (let a = 0; a < 2 * Math.PI; a += 2 * Math.PI / n) {
    particles.push(
      new Particle(
        new Vec2d(
          Math.cos(a) * r,
          Math.sin(a) * r
        )
      )
    );
  }
  return new Simulation(particles);
})();

This returns an instance of our simulated world, filled with 40 particles in a circle of radius 40.

Managing the lifetime of the simulation

We've made our world of cells, now we need to make time run. We set up a requestAnimationFrame loop (the tick(timestamp) function) that will call our simulation step to update physics as time ticks by:

// loop stuff
var lastTime = null, lastGenerated = null;

var tick = timestamp => {
  if (!lastTime) lastTime = timestamp;
  let elapsed = clamp(0, 16, timestamp - lastTime);
  lastTime = timestamp;

Here, we calculate the time that elapsed since last time tick() was called. We limit elapsed to be between 0 and 16 to keep the simulation reasonably stable.

  simulation.step(elapsed);
  update(simulation.particles);

We call step(elapsed) to advance the simulation by this number of milliseconds. We also call the update function that renders the array of particles using D3.

Making the cells divide

  if (!lastGenerated) lastGenerated = timestamp;
  // generate a new particle every 60ms
  if (simulation.particles.length < 800 && timestamp - lastGenerated > 60) {
    let idx = Math.floor(Math.random() * simulation.particles.length),
        nextIdx = (idx + 1) % simulation.particles.length,
        a = simulation.particles[idx],
        b = simulation.particles[nextIdx];
    simulation.particles.splice(
      idx,
      0,
      new Particle(
        a.position.clone().sub(b.position).scale(0.5).add(a.position)
      )
    );
    lastGenerated = timestamp;
  }

This part deals with creating a new cell within the array of particles. Every 60ms it will choose a random index within the array to insert a new cell between the chosen index and the next one. The new cell's position will be halfway between the two cells ((a + b) / 2), which is what we described above (o—o ==> o~o~o) and creates tension in the chain where the cell is inserted, because the distance between the three cells is now half of what it was before cell division.

  requestAnimationFrame(tick);
};

requestAnimationFrame(tick);

Finally, we make sure to call requestAnimationFrame when tick is done and run it for the first time.

D3 rendering

The D3 stuff is pretty straightforward for anyone used to the library, but in case you were wondering, I used the General Update Pattern that Mike Bostock recommends using.

// D3 stuff

var container = document.querySelector('#display'),
    canvas = d3.select(container)
      .append('svg')
      .attr({ width: container.offsetWidth, height: container.offsetHeight });

var update = data => {
  let w2 = container.offsetWidth / 2, h2 = container.offsetHeight / 2;
  let links = canvas.selectAll('line').data(data);
  let particles = canvas.selectAll('circle').data(data);

  // update the links
  links.enter().append('line');

  links.attr({
    x1: d => d.position.x + w2,
    y1: d => d.position.y + h2,
    x2: (d, idx) => data[(idx + 1) % data.length].position.x + w2,
    y2: (d, idx) => data[(idx + 1) % data.length].position.y + h2,
  });

  links.exit().remove();

  // update the particles
  particles.enter().append('circle');

  particles.attr({
    cx: d => d.position.x + w2,
    cy: d => d.position.y + h2,
    r: 3
  });

  particles.exit().remove();
};

Conclusion

All in all, the system described uses only a couple of very rudimentary rules, and the added randomness of cell division makes the system react in unpredictable ways. The resulting patterns, however, appear to emerge as if the "organism" were told to grow in a certain way, even though the code itself exerts no such pressure. It is all self-organisation as a consequence of the laws of nature themselves.

Your move, Creationists. :)