The tutorial codes in this page explore formation made by the mechanism inspired by cell division and growth process.

In the first code below, it defines an agent class "Cell" inheriting IParticle class to let the cell agents to move around by external or internal forces. This agent implements the following three behaviors.

• Growth : A cell agent grows its size over time. This is implemented by having a circle as the agent body geometry and increasing the radius over time in the update method.
• Repulsion : When a cell agent grows its size and starts touching other cell agents, it pushes others. This is implemented in interact method by applying force to the other cell agent if the distance of center points of two agents is smaller than summation of both radii.
• Division : Cell division is simulated by creating a new child cell by a parent cell agent and making both size half. The center positions of both cell agents are moved by the half radius towards opposite directions each other. This is implemented in update method.

The code below also control the timing of growth and division in update method by growthInterval and divisionInterval and also stop them after growthDuration time frame. The division is also controlled to happen randomly in 50 percent probability.

import igeo.*;
import processing.opengl.*;

void setup(){
  size(480,360,IG.GL);
  new Cell(new IVec(0,0,0),10).clr(1.0,0,0);
}

class Cell extends IParticle{
  int growthDuration = 2000;  //duration of growth and division
  int growthInterval = 10;
  int divisionInterval = 100;
  double radius;
  ICircle circle;
  
  Cell(IVec pos, double rad){
    super(pos, new IVec(0,0,0));
    radius = rad;
    fric(0.1);
  }
  
  void interact(ArrayList < IDynamics > agents){
    for(int i=0; i < agents.size(); i++){
      if(agents.get(i) instanceof Cell){
        Cell c = (Cell)agents.get(i);
        if(c != this){
          // push if closer than two radii
          if(c.pos().dist(pos()) < c.radius+radius){
            IVec dif = c.pos().dif(pos());
            dif.len( ((c.radius+radius)-dif.len())*100 + 50); // the closer the harder to push
            c.push(dif); 
          }
        }
      }
    }
  }
  
  void update(){
    if(IG.time() < growthDuration){
      if(time() > 0 && time()%divisionInterval==0){
        if(IRand.pct(50)){ // cell division
          IVec dir = IRand.dir(IG.zaxis);
          radius *= 0.5; //make both cell size half
          dir.len(radius);
          Cell child = new Cell(pos().cp().add(dir), radius);
          child.hsb(hue()+IRand.get(-.1,.1),saturation()+IRand.get(-.1,.1),brightness()+IRand.get(-.1,.1));
          pos().sub(dir);
        }
      }
      if(time()%growthInterval==0){
        radius += 0.1; // growing cell size
      }
    }
    // update geometry
    if(circle==null){ circle = new ICircle(pos(), radius).clr(this); }
    else{ circle.center(pos()).radius(radius); }
  }
}

The code below only reorganize the first code above by extracting part of the code of division growth from update method and putting them into grow method and divide method. The grow method also gets maxRadius parameter not to grow the size too much.

import igeo.*;
import processing.opengl.*;

void setup(){
  size(480,360,IG.GL);
  new Cell(new IVec(0,0,0), 10).clr(1.0,0,0);
}

class Cell extends IParticle{
  int growthDuration = 2000; //duration of growth and division
  int growthInterval = 10;
  int divisionInterval = 100;
  double maxRadius = 100;
  double growthSpeed = 0.1;
  
  double radius;
  ICircle circle;
  
  Cell(IVec pos, double rad){
    super(pos, new IVec(0,0,0));
    radius = rad;
    fric(0.1);
  }
  
  void interact(ArrayList < IDynamics > agents){
    for(int i=0; i < agents.size(); i++){
      if(agents.get(i) instanceof Cell){
        Cell c = (Cell)agents.get(i);
        if(c != this){
          // push if closer than two radii
          if(c.pos().dist(pos()) < c.radius+radius){
            IVec dif = c.pos().dif(pos());
            dif.len( ((c.radius+radius)-dif.len())*100+50); // the closer the harder to push
            c.push(dif);
          }
        }
      }
    }
  }
  
  void update(){
    if(IG.time() < growthDuration){
      if(time() > 0 && time()%divisionInterval==0){
        if(IRand.pct(50)){
          divide();
        }
      }
      if(time()%growthInterval==0){
        grow();
      }
    }
    // update geometry
    if(circle==null){ circle = new ICircle(pos(), radius).clr(clr()); }
    else{ circle.center(pos()).radius(radius); }
  }
  
  void divide(){ // cell division
    IVec dir = IRand.dir(IG.zaxis);
    radius *= 0.5; //make both cell size half
    dir.len(radius);
    Cell child = new Cell(pos().cp().add(dir), radius);
    child.hsb(hue()+IRand.get(-.1,.1),saturation()+IRand.get(-.1,.1),brightness()+IRand.get(-.1,.1));
    pos().sub(dir);
  }
  
  void grow(){ // growing cell size
    if(radius < maxRadius){
      radius += growthSpeed;
    }
  }
}

The next several codes show variations in controlling when cell division happens. The code below controls division by cell's location. When y coordinates of the center point of a cell is larger, it's more likely to divide because the probability is calculated by the y coordinates.

import igeo.*;
import processing.opengl.*;

void setup(){
  size(480,360,IG.GL);
  new Cell(new IVec(0,0,0), 10).clr(1.0,0,0);
}

class Cell extends IParticle{
  int growthDuration = 2000; //duration of growth and division
  int growthInterval = 10;
  int divisionInterval = 100;
  double maxRadius = 100;
  double growthSpeed = 0.1;
  
  double radius;
  ICircle circle;
  
  Cell(IVec pos, double rad){
    super(pos, new IVec(0,0,0));
    radius = rad;
    fric(0.1);
  }
  
  void interact(ArrayList < IDynamics > agents){
    for(int i=0; i < agents.size(); i++){
      if(agents.get(i) instanceof Cell){
        Cell c = (Cell)agents.get(i);
        if(c != this){
          // push if closer than two radii
          if(c.pos().dist(pos()) < c.radius+radius){
            IVec dif = c.pos().dif(pos());
            dif.len(((c.radius+radius)-dif.len())*100+50); // the closer the harder to push
            c.push(dif);
          }
        }
      }
    }
  }
  
  void update(){
    if(IG.time() < growthDuration){
      if(time() > 0 && time()%divisionInterval==0){
        double probability = (pos().y()+50)*0.8; // the upper the more probable to divide
        if(IRand.pct(probability)){
          divide();
        }
      }
      if(time()%growthInterval==0){
        grow();
      }
    }
    // update geometry
    if(circle==null){ circle = new ICircle(pos(), radius).clr(clr()); }
    else{ circle.center(pos()).radius(radius); }
  }
  
  void divide(){ // cell division
    IVec dir = IRand.dir(IG.zaxis);
    radius *= 0.5; //make both cell size half
    dir.len(radius);
    Cell child = new Cell(pos().cp().add(dir), radius);
    child.hsb(hue()+IRand.get(-.1,.1),saturation()+IRand.get(-.1,.1),brightness()+IRand.get(-.1,.1));
    pos().sub(dir);
  }
  
  void grow(){ // growing cell size
    if(radius < maxRadius){
      radius += growthSpeed;
    }
  }
}

The next code controls division by the size of a cell agent. When a cell's size (radius) becomes too big (> 5), then it divides and make the size half. It also randomize the speed of growth (growthSpeed) in each instance of cell agent.

import igeo.*;
import processing.opengl.*;

void setup(){
  size(480,360,IG.GL);
  new Cell(new IVec(0,0,0), 10);
}

class Cell extends IParticle{
  int growthDuration = 2000; //duration of growth and division
  int growthInterval = 10;
  int divisionInterval = 100;
  double maxRadius = 100;
  double growthSpeed = 0.1;
  
  double radius;
  ICircle circle;
  
  Cell(IVec pos, double rad){
    super(pos, new IVec(0,0,0));
    radius = rad;
    growthSpeed = IRand.get(0.05,0.1); 
    hsb(growthSpeed*10, 1, 1);
    fric(0.1);
  }
  
  void interact(ArrayList < IDynamics > agents){
    for(int i=0; i < agents.size(); i++){
      if(agents.get(i) instanceof Cell){
        Cell c = (Cell)agents.get(i);
        if(c != this){
          // push if closer than two radii
          if(c.pos().dist(pos()) < c.radius+radius){
            IVec dif = c.pos().dif(pos());
            dif.len(((c.radius+radius)-dif.len())*100+50); // the closer the harder to push
            c.push(dif);
          }
        }
      }
    }
  }
  
  void update(){
    if(IG.time() < growthDuration){
      if(radius > 5){ // divide when bigger than 5
        divide();
      }
      if(time()%growthInterval==0){
        grow();
      }
    }
    // update geometry
    if(circle==null){ circle = new ICircle(pos(), radius).clr(clr()); }
    else{ circle.center(pos()).radius(radius); }
  }
  
  void divide(){ // cell division
    IVec dir = IRand.dir(IG.zaxis);
    radius *= 0.5; //make both cell size half
    dir.len(radius);
    Cell child = new Cell(pos().cp().add(dir), radius);
    pos().sub(dir);
  }
  
  void grow(){ // growing cell size
    if(radius < maxRadius){
      radius += growthSpeed;
    }
  }
}

The next example shows a control of division by attractor points. With the analogy of food or nutrition gives cells energy to grow and divide, the attractor is implemented as Nutrition class in the code and when cells location is closer than a certain distance (Nutrition's radius), the cells divide themselves. With another analogy of nutrition being consumed, the size of a Nutrition agent decreases when it lets a cell agent divide.

Cell agent class also adds active boolean field to control division any time before update method and Nutrition agent turn cell's active field to true to activate division in its feed method. active field is set to be false in divide method to reset the activation every time.

import igeo.*;
import processing.opengl.*;

void setup(){
  size(480,360,IG.GL);
  for(int i=0; i < 8; i++){
    new Nutrition(IRand.pt(-100,-100,100,100), 20);
  }
  for(int i=0; i < 20; i++){
    new Cell(IRand.pt(-100,-100,100,100), 10);
  }
}

class Nutrition extends IAgent{
  IVec pos;
  double radius;
  ICircle circle;
  Nutrition(IVec p, double rad){
    pos = p;
    radius = rad;
    circle = new ICircle(pos, radius).clr(1.,0,0).weight(2);
  }
  
  void feed(Cell c){
    if(pos.dist(c.pos()) < radius){
      c.active = true; // activate cell
      radius -= 0.05; // shrink 
      circle.radius(radius);
    }
  } 
}

class Cell extends IParticle{
  int growthDuration = 1400; //duration of growth and division
  int growthInterval = 10;
  int divisionInterval = 100;
  double maxRadius = 100;
  double growthSpeed = 0.1;
  
  double radius;
  ICircle circle;
  boolean active = false;
  
  Cell(IVec pos, double rad){
    super(pos, new IVec(0,0,0));
    radius = rad;
    fric(0.1);
  }
  
  void interact(ArrayList < IDynamics > agents){
    for(int i=0; i < agents.size(); i++){
      if(agents.get(i) instanceof Cell){
        Cell c = (Cell)agents.get(i);
        if(c != this){
          // push if closer than two radii
          if(c.pos().dist(pos()) < c.radius+radius){
            IVec dif = c.pos().dif(pos());
            dif.len(((c.radius+radius)-dif.len())*100+50); // the closer the harder to push
            c.push(dif);
          }
        }
      }
      if(agents.get(i) instanceof Nutrition){
        if(time()%divisionInterval==0 && IG.time() < growthDuration){
          Nutrition n = (Nutrition)agents.get(i);
          n.feed(this);
        }
      }
    }
  }
  
  void update(){
    if(IG.time() < growthDuration){
      if(time()>0&&time()%divisionInterval==0 ){
        if(active){ // divide when active flag is on
          divide();
        }
      }
      if(time()%growthInterval==0){
        grow();
      }
    }
    // update geometry
    hsb(radius/10,1,0.8); 
    if(circle==null){ circle = new ICircle(pos(), radius).clr(clr()); }
    else{ circle.center(pos()).radius(radius).clr(clr()); }
  }
  
  void divide(){ // cell division
    IVec dir = IRand.dir(IG.zaxis);
    radius *= 0.5; //make both cell size half
    dir.len(radius);
    Cell child = new Cell(pos().cp().add(dir), radius);
    pos().sub(dir);
    active = false; //reset after division
  }
  
  void grow(){ // growing cell size
    if(radius < maxRadius){
      radius += growthSpeed;
    }
  }  
}

The next two codes show an example of division control by measuring how crowded it is around a cell. The cell agent measures how crowded around it by counting number of other cell agents within a certain distance range (neighborDist) in interact method. If the total number of neighbors (neighborCount) within the distance range is not too many (<= 3), the cell activates division by turning active variable true in the code below.

import igeo.*;
import processing.opengl.*;

void setup(){
  size(480,360,IG.GL);
  new Cell(new IVec(0,0,0), 1).clr(1.0,0,0);
}

class Cell extends IParticle{
  int growthDuration = 1800; //duration of growth and division
  int growthInterval = 10;
  int divisionInterval = 100;
  double maxRadius = 100;
  double growthSpeed = 0.1;
  
  double radius;
  ICircle circle;
  boolean active = false;
  
  Cell(IVec pos, double rad){
    super(pos, new IVec(0,0,0));
    radius = rad;
    fric(0.1);
  }
  
  void interact(ArrayList < IDynamics > agents){
    int neighborCount=0;
    double neighborDist = 1;
    for(int i=0; i < agents.size(); i++){
      if(agents.get(i) instanceof Cell){
        Cell c = (Cell)agents.get(i);
        if(c != this){
          // push if closer than two radii
          if(c.pos().dist(pos()) < c.radius+radius){
            IVec dif = c.pos().dif(pos());
            dif.len(((c.radius+radius)-dif.len())*100+50); // the closer the harder to push
            c.push(dif);
          }
          if(c.pos().dist(pos()) < c.radius+radius+neighborDist){
            neighborCount++; // count close neighbors
          }
        }
      }
    }
    // activate when not many neighbors
    if(neighborCount <= 3){
      active = true;
    }
    else{
      active = false;
    }
  }
  
  void update(){
    if(IG.time() < growthDuration){
      if(time() > 0 && time()%divisionInterval==0){
        if(active){ // divide when active flag is on
          divide();
        }
      }
      if(time()%growthInterval==0){
        grow();
      }
    }
    // update geometry
    if(circle==null){ circle = new ICircle(pos(), radius).clr(clr()); }
    else{ circle.center(pos()).radius(radius); }
  }
  
  void divide(){ // cell division
    IVec dir = IRand.dir(IG.zaxis);
    radius *= 0.5; //make both cell size half
    dir.len(radius);
    Cell child = new Cell(pos().cp().add(dir), radius);
    child.hsb(hue()+IRand.get(-.1,.1),saturation()+IRand.get(-.1,.1),brightness()+IRand.get(-.1,.1));
    pos().sub(dir);
    active=false; //reset after division
  }
  
  void grow(){ // growing cell size
    if(radius < maxRadius){
      radius += growthSpeed;
    }
  }
}

The code below has an opposite rule of the previous one in the division control. Instead of activating division when a cell sees small number of neighbors, it activates division when a cell sees many neighbors (> 5). Because it starts with one cell, cells are randomly activated for division in the first 600 time frame, and then after that, it activates only when it finds more than 5 neighbors.

import igeo.*;
import processing.opengl.*;

void setup(){
  size(480,360,IG.GL);
  new Cell(new IVec(0,0,0), 10).clr(1.0,0,0);
}

class Cell extends IParticle{
  int growthDuration = 1000; //duration of growth and division
  int growthInterval = 10;
  int divisionInterval = 50;
  double maxRadius = 100;
  
  double growthSpeed = 0.2;
  double radius;
  ICircle circle;
  boolean active = false;
  
  Cell(IVec pos, double rad){
    super(pos, new IVec(0,0,0));
    radius = rad;
    fric(0.1);
  }
  
  void interact(ArrayList < IDynamics > agents){
    int neighborCount=0;
    double neighborDist = 2; 
    for(int i=0; i < agents.size(); i++){
      if(agents.get(i) instanceof Cell){
        Cell c = (Cell)agents.get(i);
        if(c != this){
          // push if closer than two radii
          if(c.pos().dist(pos()) < c.radius+radius){
            IVec dif = c.pos().dif(pos());
            dif.len(((c.radius+radius)-dif.len())*100+50); // the closer the harder to push
            c.push(dif);
          }
          if(c.pos().dist(pos()) < c.radius+radius+neighborDist){
            neighborCount++; // count close neighbors
          }
        }
      }
    }
    if(IG.time() >= 600){
      // activate when it has many neighbors
      if(neighborCount > 5){
        active = true;
      }
      else{
        active = false;
      }
    }
  }
  
  void update(){
    if(IG.time() < growthDuration){
      if(time() > 0 && time()%divisionInterval==0 ){
        if(IG.time() < 600){ // randomly activate in early stage
          if(IRand.pct(50)){
            active = true; 
          }
        }
        if(active){ // divide when active flag is on
          divide();
        }
      }
      if(time()%growthInterval==0){
        grow();
      }
    }
    // update geometry
    if(circle==null){ circle = new ICircle(pos(), radius).clr(clr()); }
    else{ circle.center(pos()).radius(radius); }
  }
  
  void divide(){ // cell division
    IVec dir = IRand.dir(IG.zaxis);
    radius *= 0.5; //make both cell size half
    dir.len(radius);
    Cell child = new Cell(pos().cp().add(dir), radius);
    child.hsb(hue()+IRand.get(-.1,.1),saturation()+IRand.get(-.1,.1),brightness()+IRand.get(-.1,.1));
    pos().sub(dir);
    active = false; //reset after division
  }
  
  void grow(){ // growing cell size
    if(radius < maxRadius){
      radius += growthSpeed;
    }
  }
}

Many cells in nature have a mechanism to adhere to other cells. This cell adhesion or cell junction mechanism is simulated by connecting two cells with spring-like tensile force. This connection is created when a cell divides itself into two. The parent cell put a spring link between itself and a new divided child.

The code below introduces a new class CellLink to implement the connection and adhesion. To simulate the adhesion by spring-like force, CellLink has an interact method to calculate a force vector from a difference vector of two positions of the linked cells and pull them if they are farther than the summation of two radii or push them if they are too close. An instance of CellLink is created at Cell class's divide method by passing the cell itself and a new child cell.

import igeo.*;
import processing.opengl.*;

void setup(){
  size(480,360,IG.GL);
  new Cell(new IVec(0,0,0), 10).clr(0,0,1.0);
}

class Cell extends IParticle{
  int growthDuration = 1800; //duration of growth and division
  int growthInterval = 10;
  int divisionInterval = 100;
  double maxRadius = 100;
  double growthSpeed = 0.1;
  
  double radius;
  ICircle circle;
  boolean active=false;
  
  Cell(IVec pos, double rad){
    super(pos, new IVec(0,0,0));
    radius = rad;
    fric(0.1);
  }
  
  void interact(ArrayList< IDynamics > agents){
    for(int i=0; i < agents.size(); i++){
      if(agents.get(i) instanceof Cell){
        Cell c = (Cell)agents.get(i);
        if(c != this){
          // push if closer than two radii
          if(c.pos().dist(pos()) < c.radius+radius){
            IVec dif = c.pos().dif(pos());
            dif.len(((c.radius+radius)-dif.len())*100+50); // the closer the harder to push
            c.push(dif);
          }
        }
      }
    }
  }
  
  void update(){
    if(IG.time() < growthDuration){
      if(time()>0 && time()%divisionInterval==0){
        if(IRand.pct(50)){ // random division
          active = true;
        }
        if(active){ // divide when active flag is on
          divide();
        }
      }
      if(time()%growthInterval==0){
        grow();
      }
    }
    // update geometry
    if(circle==null){ circle = new ICircle(pos(), radius).clr(clr()); }
    else{ circle.center(pos()).radius(radius); }
  }
  
  void divide(){ // cell division
    IVec dir = IRand.dir(IG.zaxis);
    radius *= 0.5; //make both cell size half
    dir.len(radius);
    Cell child = new Cell(pos().cp().add(dir), radius);
    child.hsb(hue()+IRand.get(-.1,.1),saturation()+IRand.get(-.1,.1),brightness()+IRand.get(-.1,.1));
    pos().sub(dir);
    active = false; //reset after division
    new CellLink(this, child);
  }
  
  void grow(){ // growing cell size
    if(radius < maxRadius){
      radius += growthSpeed;
    }
  }
}

class CellLink extends IAgent{
  Cell cell1, cell2;
  ICurve line;
  CellLink(Cell c1, Cell c2){
    cell1 = c1; cell2 = c2;
    line = new ICurve(c1.pos(), c2.pos()).clr(1.0,0,0);
  }
  
  void interact(ArrayList< IDynamics > agents){
    // spring force
    IVec dif = cell1.pos().dif(cell2.pos());
    double force = (dif.len()-(cell1.radius+cell2.radius))/(cell1.radius+cell2.radius)*200;
    dif.len(force);
    cell1.pull(dif);
    cell2.push(dif);
  }
}

Just to show the effect of tensile links and the network structure of links clearer, the code below adds another repulsion behavior in interact method to push other cells when they are within a certain distance range (30) in addition to the original repulsion force not to have their body defined by radius overlapped.

import igeo.*;
import processing.opengl.*;

void setup(){
  size(480,360,IG.GL);
  new Cell(new IVec(0,0,0), 10).clr(0,0,1.0);
}

class Cell extends IParticle{
  int growthDuration = 1800; //duration of growth and division
  int growthInterval = 10;
  int divisionInterval = 100;
  double maxRadius = 100;
  double growthSpeed = 0.1;
  
  double radius;
  ICircle circle;
  boolean active = false;
  
  Cell(IVec pos, double rad){
    super(pos, new IVec(0,0,0));
    radius = rad;
    fric(0.1);
  }
  
  void interact(ArrayList< IDynamics > agents){
    for(int i=0; i < agents.size(); i++){
      if(agents.get(i) instanceof Cell){
        Cell c = (Cell)agents.get(i);
        if(c != this){
          // push if closer than two radii
          if(c.pos().dist(pos()) < c.radius+radius){
            IVec dif = c.pos().dif(pos());
            dif.len(((c.radius+radius)-dif.len())*100+50); // the closer the harder to push
            c.push(dif);
          }
          // additional repulsion for smoother organization
          if(c.pos().dist(pos()) < c.radius+radius+30){
            IVec dif = c.pos().dif(pos());
            dif.len(1); // constant force
            c.push(dif);
          }
        }
      }
    }
  }
  
  void update(){
    if(IG.time() < growthDuration){
      if(time()>0 && time()%divisionInterval==0){
        if(IRand.pct(50)){ // random division
          active = true;
        } 
        if(active){ // divide when active flag is on
          divide();
        }
      }
      if(time()%growthInterval==0){
        grow();
      }
    }
    // update geometry
    if(circle==null){ circle = new ICircle(pos(), radius).clr(clr()); }
    else{ circle.center(pos()).radius(radius); }
  }
  
  void divide(){ // cell division
    IVec dir = IRand.dir(IG.zaxis);
    radius *= 0.5; //make both cell size half
    dir.len(radius);
    Cell child = new Cell(pos().cp().add(dir), radius);
    child.hsb(hue()+IRand.get(-.1,.1),saturation()+IRand.get(-.1,.1),brightness()+IRand.get(-.1,.1));
    pos().sub(dir);
    active = false; //reset after division
    new CellLink(this, child);
  }
  
  void grow(){ // growing cell size
    if(radius < maxRadius){
      radius += growthSpeed;
    }
  }
}

class CellLink extends IAgent{
  Cell cell1, cell2;
  ICurve line;
  CellLink(Cell c1, Cell c2){
    cell1 = c1; cell2 = c2;
    line = new ICurve(c1.pos(), c2.pos()).clr(1.,0,0);
  }
  
  void interact(ArrayList< IDynamics > agents){
    // spring force
    IVec dif = cell1.pos().dif(cell2.pos());
    double force = (dif.len()-(cell1.radius+cell2.radius))/(cell1.radius+cell2.radius)*200;
    dif.len(force);
    cell1.pull(dif);
    cell2.push(dif);
  }
}

Instead of random activation of division or activation depending on surrounding condition at the moment, the next two codes show a division control depending on the history of parent and child through time. In the cascading activation, a parent cell activates a child cell but the parent cell deactivates itself. When the first cell is activated in setup method, this activation cascades to the youngest child cell and there is always only one active cell at a moment.

import igeo.*;
import processing.opengl.*;

void setup(){
  size(480,360,IG.GL);
  Cell cell = new Cell(new IVec(0,0,0),10);
  cell.clr(0,0,1.0);
  cell.active = true;
}

class Cell extends IParticle{
  int growthDuration = 2000; //duration of growth and division
  int growthInterval = 10;
  int divisionInterval = 50;
  double maxRadius = 100;
  double growthSpeed = 0.1;
  
  double radius;
  ICircle circle;
  boolean active = false;
  
  Cell(IVec pos, double rad){
    super(pos, new IVec(0,0,0));
    radius = rad;
    fric(0.1);
  }
  
  void interact(ArrayList< IDynamics > agents){
    for(int i=0; i < agents.size(); i++){
      if(agents.get(i) instanceof Cell){
        Cell c = (Cell)agents.get(i);
        if(c != this){
          // push if closer than two radii
          if(c.pos().dist(pos()) < c.radius+radius){
            IVec dif = c.pos().dif(pos());
            dif.len(((c.radius+radius)-dif.len())*100+50); // the closer the harder to push
            c.push(dif);
          }
          // additional repulsion for smoother organization
          if(c.pos().dist(pos()) < c.radius+radius+100){
            IVec dif = c.pos().dif(pos());
            dif.len(1); // constant force
            c.push(dif);
          }
        }
      }
    }
  }
  
  void update(){
    if(IG.time() < growthDuration){
      if(time() > 0 && time()%divisionInterval==0){
        if(active){ // divide when active flag is on
          divide();
        }
      }
      if(time()%growthInterval==0){
        grow();
      }
    }
    // update geometry
    if(circle==null){ circle = new ICircle(pos(), radius).clr(clr()); }
    else{ circle.center(pos()).radius(radius); }
  }
  
  void divide(){ // cell division
    IVec dir = IRand.dir(IG.zaxis);
    radius *= 0.5; //make both cell size half
    dir.len(radius);
    Cell child = new Cell(pos().cp().add(dir), radius);
    child.hsb(hue()+IRand.get(-.1,.1),saturation()+IRand.get(-.1,.1),brightness()+IRand.get(-.1,.1));
    pos().sub(dir);
    child.active = true; //activate child
    active = false; //deactivate itself
    new CellLink(this, child);
  }
  
  void grow(){ // growing cell size
    if(radius < maxRadius){
      radius += growthSpeed;
    }
  }
}

class CellLink extends IAgent{
  Cell cell1, cell2;
  ICurve line;
  CellLink(Cell c1, Cell c2){
    cell1 = c1; cell2 = c2;
    line = new ICurve(c1.pos(), c2.pos()).clr(1.0,0,0);
  }
  
  void interact(ArrayList< IDynamics > agents){
    // spring force
    IVec dif = cell1.pos().dif(cell2.pos());
    double force = (dif.len()-(cell1.radius+cell2.radius))/(cell1.radius+cell2.radius)*200;
    dif.len(force);
    cell1.pull(dif);
    cell2.push(dif);
  }
}

The code below shows another parent-child history dependent activation control of cell division. Instead of a parent cell deactivate itself, it maintains activation after creating a new child and activating the child. Because the still active parent cell will create another child cell in the next time, the parent cell creates multiple children and branched links. Because if cells create child cells in every cycle, they will produce exponentially large number of cells, the code below limits only 3% of parent cells who have two links (back and front) to maintain the activation.

import igeo.*;
import processing.opengl.*;

void setup(){
  size(480,360,IG.GL);
  Cell cell = new Cell(new IVec(0,0,0), 10);
  cell.clr(0,0,1.0);
  cell.active = true;
}

class Cell extends IParticle{
  int growthDuration = 1000; //duration of growth and division
  int growthInterval = 10;
  int divisionInterval = 50;
  double maxRadius = 100;
  double growthSpeed = 0.1;
  
  double radius;
  ICircle circle;
  boolean active = false;
  
  Cell(IVec pos, double rad){
    super(pos, new IVec(0,0,0));
    radius = rad;
    fric(0.1);
  }
  
  void interact(ArrayList< IDynamics > agents){
    for(int i=0; i < agents.size(); i++){
      if(agents.get(i) instanceof Cell){
        Cell c = (Cell)agents.get(i);
        if(c != this){
          // push if closer than two radii
          if(c.pos().dist(pos()) < c.radius+radius){
            IVec dif = c.pos().dif(pos());
            dif.len(((c.radius+radius)-dif.len())*100+50); // the closer the harder to push
            c.push(dif);
          }
          // additional repulsion for smoother organization
          if(c.pos().dist(pos()) < c.radius+radius+50){
            IVec dif = c.pos().dif(pos());
            dif.len(1); // constant force
            c.push(dif);
          }
        }
      }
    }
  }
  
  void update(){
    if(IG.time() < growthDuration){
      if(time()>0 && time()%divisionInterval==0){
        if(active){ // divide when active flag is on
          divide();
        }
      }
      if(time()%growthInterval==0){
        grow();
      }
    }
    // update geometry
    if(circle==null){ circle = new ICircle(pos(), radius).clr(clr()); }
    else{ circle.center(pos()).radius(radius); }
  }
  
  void divide(){ // cell division
    IVec dir = IRand.dir(IG.zaxis);
    radius *= 0.5; //make both cell size half
    dir.len(radius);
    Cell child = new Cell(pos().cp().add(dir), radius);
    child.hsb(hue()+IRand.get(-.1,.1),saturation()+IRand.get(-.1,.1),brightness()+IRand.get(-.1,.1));
    pos().sub(dir);
    child.active =true; //activate child
    if(IRand.pct(70)){ //some deactivated others stay active
      active = false; //deactivate itself
    }
    new CellLink(this, child);
  }
  
  void grow(){ // growing cell size
    if(radius < maxRadius){
      radius += growthSpeed;
    }
  }
  
}

class CellLink extends IAgent{
  Cell cell1, cell2;
  ICurve line;
  CellLink(Cell c1, Cell c2){
    cell1 = c1; cell2 = c2;
    line = new ICurve(c1.pos(), c2.pos()).clr(0,0.5,0);
  }
  
  void interact(ArrayList< IDynamics > agents){
    // spring force
    IVec dif = cell1.pos().dif(cell2.pos());
    double force = (dif.len()-(cell1.radius+cell2.radius))/(cell1.radius+cell2.radius)*200;
    dif.len(force);
    cell1.pull(dif);
    cell2.push(dif);
  }
}

The code below produces branching link network similar to the previous one but its division control algorithm is different. A cell is activated for division only when it has two links (front and back) and in a probability of 3%. To know the number of links connected to the cell, the cell needs to remember what links it has got so far. To keep track of CellLink instances, Cell class adds a new field links which is a variable-length array to contain all CellLink instances connected to the cell and a CellLink instance adds itself to the list in both sides of cells in the constructor CellLink.

import igeo.*;
import processing.opengl.*;

void setup(){
  size(480,360,IG.GL);
  Cell cell = new Cell(new IVec(0,0,0), 10);
  cell.clr(0,0,1.0);
  cell.active = true;
}

class Cell extends IParticle{
  int growthDuration = 1800; //duration of growth and division
  int growthInterval = 10;
  int divisionInterval = 50;
  double maxRadius = 100;
  double growthSpeed = 0.1;
  
  ArrayList< CellLink > links; //store links
  double radius;
  ICircle circle;
  boolean active = false;
  
  Cell(IVec pos, double rad){
    super(pos, new IVec(0,0,0));
    radius = rad;
    links = new ArrayList< CellLink >();
    fric(0.1);
  }
  
  void interact(ArrayList< IDynamics > agents){
    for(int i=0; i < agents.size(); i++){
      if(agents.get(i) instanceof Cell){
        Cell c = (Cell)agents.get(i);
        if(c != this){
          // push if closer than two radii
          if(c.pos().dist(pos()) < c.radius+radius){
            IVec dif = c.pos().dif(pos());
            dif.len(((c.radius+radius)-dif.len())*100+50); // the closer the harder to push
            c.push(dif);
          }
          // additional repulsion for smoother organization
          if(c.pos().dist(pos()) < c.radius+radius+100){
            IVec dif = c.pos().dif(pos());
            dif.len(1); // constant force
            c.push(dif);
          }
        }
      }
    }
  }
  
  void update(){
    if(IG.time() < growthDuration){
      if(time()>0 && time()%divisionInterval==0){
        if(links.size()==2 && IRand.pct(3)){
          active = true;
        }
        if(active){ // divide when active flag is on
          divide();
        }
      }
      if(time()%growthInterval==0){
        grow();
      }
    }
    // update geometry
    if(circle==null){ circle = new ICircle(pos(), radius).clr(clr()); }
    else{ circle.center(pos()).radius(radius); }
  }
  
  void divide(){ // cell division
    IVec dir = IRand.dir(IG.zaxis);
    radius *= 0.5; //make both cell size half
    dir.len(radius);
    Cell child = new Cell(pos().cp().add(dir), radius);
    child.hsb(hue()+IRand.get(-.1,.1),saturation()+IRand.get(-.1,.1),brightness()+IRand.get(-.1,.1));
    pos().sub(dir);
    child.active =true; //activate child
    active = false; //deactivate itself
    new CellLink(this, child);
  }
  
  void grow(){ // growing cell size
    if(radius < maxRadius){
      radius += growthSpeed;
    }
  }
  
}

class CellLink extends IAgent{
  Cell cell1, cell2;
  ICurve line;
  CellLink(Cell c1, Cell c2){
    cell1 = c1; cell2 = c2;
    cell1.links.add(this); // register this link to cells
    cell2.links.add(this);
    line = new ICurve(c1.pos(), c2.pos()).clr(0,0.5,0);
  }
  
  void interact(ArrayList< IDynamics > agents){
    // spring force
    IVec dif = cell1.pos().dif(cell2.pos());
    double force = (dif.len()-(cell1.radius+cell2.radius))/(cell1.radius+cell2.radius)*200;
    dif.len(force);
    cell1.pull(dif);
    cell2.push(dif);
  }
}

The tutorial code shows a different way to add a new link to a child cell and it creates a closed loop of links. So far a parent cell just put one link between a child cell and the parent and the newborn child cell always connected to only one cell. The diagram below shows another way to connect. If a parent cell already has one or more links, a new child cell can be inserted in the middle of one existing link.

When a child cell is inserted at an existing link, the code actually deletes the existing link and create a new links from the parent to the child and another link from the child to the cell at the opposite end of the existing link. To delete a link, CellLink class adds a method delete and this method unregisters the link from the list of links inside each cell. For a cell to find the cell at the opposite end of link, CellLink class also adds a method oppositeCell.

In the method divide in Cell class, the part of the code to create a child cell is extracted and put inside a new method createChild. divide method adds more codes to change the behavior depending on the number of links in the cell. When a cell has no link, it just creates one link to the new child which is the same behavior with previous codes. If the cell has one link, which means this cell is at open end of links, it creates two links from the cell to its child cell, and also from the child to the opposite end of the existing link. This creates a closed triangular links. If the number of links is two, then it performs the insertion behavior by first selecting which of two links to insert and second deleting the link and the third put two links from the parent to the child and the child to the opposite of the deleted link.

The code below also adds another behavior in interact method. It's a type of repulsion behavior which is same with the separation behavior in Boid algorithm. The behavior lets the cell to move away from the center of neighbors.

import igeo.*;
import processing.opengl.*;

void setup(){
  size(480,360,IG.GL);
  IConfig.syncDrawAndDynamics=true; //not to crash when some geometry is deleted while drawing
  new Cell(new IVec(0,0,0), 10).clr(0,0,1.0);
}

class Cell extends IParticle{
  int growthDuration = 1000; //duration of growth and division
  int growthInterval = 10;
  int divisionInterval = 50;
  double maxRadius = 100;
  double growthSpeed = 0.1;
  
  ArrayList< CellLink > links; //store links
  double radius;
  ICircle circle;
  boolean active = false;
  
  Cell(IVec pos, double rad){
    super(pos, new IVec(0,0,0));
    radius = rad;
    links = new ArrayList< CellLink >();
    fric(0.1);
  }
  
  void interact(ArrayList< IDynamics > agents){
    IVec neighborCenter = new IVec(0,0,0);
    int neighborCount=0;
    for(int i=0; i < agents.size(); i++){
      if(agents.get(i) instanceof Cell){
        Cell c = (Cell)agents.get(i);
        if(c != this){
          // push if closer than two radii
          if(c.pos().dist(pos()) < c.radius+radius){
            IVec dif = c.pos().dif(pos());
            dif.len(((c.radius+radius)-dif.len())*100+50); // the closer the harder to push
            c.push(dif);
          }
          // count neighbors and calculate their center
          if(c.pos().dist(pos()) < c.radius+radius + radius*2){
            neighborCenter.add(c.pos());
            neighborCount++;
          }
        }
      }
    }
    
    if(neighborCount > 0){ // push from center of neighbors
      neighborCenter.div(neighborCount);
      IVec dif = pos().dif(neighborCenter).len(20); // constant force
      push(dif);
    }
  }
  
  void update(){
    if(IG.time() < growthDuration){
      if(time() > 0 && time()%divisionInterval==0){
        if(IRand.pct(50)){ // random division
          active = true;
        }
        if(active){ // divide when active flag is on
          divide();
        }
      }
      if(time()%growthInterval==0){
        grow();
      }
    }
    // update geometry
    if(circle==null){ circle = new ICircle(pos(), radius).clr(this); }
    else{ circle.center(pos()).radius(radius); }
  }
  
  void divide(){ // cell division
    if(links.size()==0){ // dot state
      Cell child = createChild(IRand.dir(IG.zaxis));
      new CellLink(this, child);
    }
    else if(links.size()==1){ // line state
      Cell child = createChild(IRand.dir(IG.zaxis));
      new CellLink(child, links.get(0).cell1); // making a triangle loop
      new CellLink(child, links.get(0).cell2);
    }
    else if(links.size()==2){ // string state
      CellLink dividingLink = links.get(IRand.getInt(0,1)); // pick one link out of two
      Cell c = dividingLink.oppositeCell(this); // other cell on the link
      IVec dir = c.pos().dif(pos()); // dividing direction is link direction
      Cell child = createChild(dir);
      dividingLink.del(); // delete picked link
      new CellLink(this, child); // create two new links
      new CellLink(child, c);
    }
  }
  
  void grow(){ // growing cell size
    if(radius < maxRadius){
      radius += growthSpeed;
    }
  }
  
  Cell createChild(IVec dir){
    radius *= 0.5; //make both cell size half
    dir.len(radius);
    Cell child = new Cell(pos().cp().add(dir), radius);
    child.hsb(hue()+IRand.get(-.1,.1),saturation()+IRand.get(-.1,.1),brightness()+IRand.get(-.1,.1));
    pos().sub(dir);
    active = false; //reset activation
    return child;
  }
}

class CellLink extends IAgent{
  Cell cell1, cell2;
  ICurve line;
  CellLink(Cell c1, Cell c2){
    cell1 = c1; cell2 = c2;
    cell1.links.add(this); // register this link to cells{
    cell2.links.add(this);
    line = new ICurve(c1.pos(), c2.pos()).clr(1.0,0,0);
  }
  
  void interact(ArrayList< IDynamics > agents){
    // spring force
    IVec dif = cell1.pos().dif(cell2.pos());
    double force = (dif.len()-(cell1.radius+cell2.radius))/(cell1.radius+cell2.radius)*200;
    dif.len(force);
    cell1.pull(dif);
    cell2.push(dif);
  }
  
  void del(){
    cell1.links.remove(this); // unregister from cells
    cell2.links.remove(this);
    line.del(); // delete line geometry
    super.del(); // stop agent
  }
  
  Cell oppositeCell(Cell c){ // find other cell on the link
    if(cell1==c) return cell2;
    if(cell2==c) return cell1;
    IG.err("Link does not contain the input cell");
    return null; 
  }
}