Cell Division and Growth Algorithm: Adding springs

I think I’ve managed to compute the repulsion force but for some reason cells start to shake and overlap when they are attached to springs.

No springs, no overlapping -->

Springs: overlapping -->

This is all the more confusing that the springs rest lengths are updated at each iteration (mother radius + child radius).

I reallly need your help !

add_library('verletphysics')
add_library('toxiclibscore')
from toxi.physics2d import VerletPhysics2D, VerletParticle2D
from toxi.physics2d.behaviors import GravityBehavior

id = 0

def setup():
    global collection, physics, d
    size(690, 460, P2D)
    smooth(8)
        
    physics = VerletPhysics2D()
    physics.setWorldBounds(Rect(0, 0, width, height))
    
    
    collection = [Cell(width>>1, height>>1)]
    
    #Adding first particle to array list
    physics.addParticle(collection[0])
    
    #Setting a fixed repulsion force for this Cell only
    physics.addBehavior(AttractionBehavior(collection[0], collection[0].radius*5, -2.4))
    
    #Dictionnary to keep trace of 'mother' / 'child' relations
    d = {}
    

        
def draw():
    global physics, collection
    background(30)
    
    #Updating physics engine
    physics.update()
    
    #Locking first Cell
    collection[0].lock()
    
    
    #Making cells grow and repel each other
    for c in collection:
        c.growth()
        c.display()
        c.interact()
        
        
    #Updating rest lenghts
    if len(collection) > 1:
        for i in range(len(physics.springs)):
            physics.springs.get(i).setRestLength(collection[-2].radius + collection[-1].radius + 4)
            physics.springs.get(i).setStrength(.5)
   
            
 
    #Drawing connections
    for e in d:
        for id in d[e]:
            strokeWeight(2 - id*.006)
            stroke(80 + id*.6, 140 + id, 255 - id)
            line(collection[e].x(), collection[e].y(), collection[id].x(), collection[id].y())


    
class Cell(VerletParticle2D):
    growthInterval = 5
    divisionInterval = 50
    
    def __init__(self, x, y):
        super(Cell, self).__init__(x, y)
        self.position = Vec2D(x, y)
        self.radius = 5
        self.mass = self.radius
        self.time = 0
        self.id = 0        
        
                                
    def growth(self):
        global id
        self.time += 1
        if len(collection) < 40 :
            if self.time > 0 and self.time%Cell.divisionInterval == 0:
                
                if random(1) > .5:
                    
                    #Divide radius by 2
                    self.radius *= .5
                    
                    #Update id
                    id += 1
                    
                    #Create random vector
                    dir = Vec2D.randomVector()
            
                    #Adding random vector to child position
                    child = Cell(self.x() + dir.x(), self.y() + dir.y())
                    
                    #Setting child's id
                    child.id = id
                    
                    #Keep trace of 'mother' <-> 'child' relations
                    if self.id not in d:
                        d[self.id] = [child.id]
                    else:
                        d[self.id].append(child.id)
                    
                    #Append cell to array list
                    collection.append(child)
                    
                    #Adding cell to physics engine
                    physics.addParticle(child)
                    
                    #Adding spring between mother and child
                    s = VerletSpring2D(self, child, self.radius + child.radius, .5)
                    physics.addSpring(s)
                    
                
            
            if self.time%Cell.growthInterval == 0:
                self.radius += .5

                
    def interact(self):
        
        #Check every other cells
        for a in collection:
            if a is not self:
                
                #Distance between other cells
                d = dist(a.x(), a.y(), self.x(), self.y())
                if d < (a.radius + self.radius):
                    
                    ### Push other cells outside its perimeter ###
                                        
                    #Direction
                    dir = PVector.sub(PVector(self.x(), self.y()), PVector(a.x(), a.y()))
                    dis = dir.magSq()
                    dis = constrain(dis, 10, (self.radius + a.radius))
                    dir.normalize()
                    
                    #Magnitude
                    m = (self.mass * a.mass) / dis
                    
                    #Combining magnitude and direction
                    dir.mult(-m)
                    
                    #Converting to Vec2D + adding force to physics engine
                    dir = Vec2D(dir.x, dir.y)
                    a.addForce(dir)
                   
                                  
    def display(self):
        stroke(255)
        strokeWeight(3)
        point(self.x(), self.y())
        fill(300 - self.id*.4, 80 + self.id*.3, 120 + self.id, 10)
        strokeWeight(.7)
        stroke(300 - self.id*.5, 80 + self.id*.3, 120 + self.id, map(self.time, 0, 20, 0, 255))
        ellipse(self.x(), self.y(), self.radius*2, self.radius*2)