Model division of cancer cells on a grid

Question

I have a 5000x5000 grid, and I'm trying to implement a simple model of cancer division in MATLAB. Initially, it picks a random point (x,y) and makes that cell a cancer cell. On the first iteration, it divides - the parent cell stays in it's place, the daughter cell is randomly assigned to any neighbouring cell.
Easy so far.

My problem is this: on successive iterations, a daughter cell will often be assigned to a cell that already has a cancer cell. In this case, I want the daughter cell to take its place and "bump" the cell already there to an adjacent cell. If that adjacent cell is empty, it is filled and the process stops. If not, the cell already in that place is bumped and so on until the last cell finds an empty space and the process stops.

This should be simple, but I have no idea how to code it up and what kind of loops to use.
I'm a physical scientists rather than a programmer, so please treat me like a simpleton!

Answer 1

Here is a function I hacked together that roughly meets the specs you provided.

I does slow down as the number of cancerous cells gets large.

Basically I have a few variables, the NxN matrix that represents the grid of cell locations (i call this a plate as grid is the name of an existing matlab function)

A vector of points that I can iterate through quickly. I pick a seed location and then run a while loop until the grid is full.

On each loop iteration I perform the following for each cell:

• Generate a random number to determine if that cell should divide
• Generate a random direction to divide
• Find the first open plate position in that direction
• Populate that position

I haven't tested it extensively but it appears to work.

function simulateCancer(plateSize, pDivide)

plate = zeros(plateSize, plateSize);
nCells = 1;
cellLocations = zeros(plateSize*plateSize,2);

initX = randi(plateSize);
initY = randi(plateSize);

cellLocations(nCells,:) = [initX, initY];

plate(initX, initY) = 1;

f = figure;
a = axes('Parent', f);
im = imagesc(plate, 'Parent', a);


while(nCells < (plateSize * plateSize))
    currentGeneration = currentGeneration+1;
    for i = 1:nCells
        divide = rand();
        if divide <= pDivide
            divideLocation = cellLocations(i,:);
            divideDir = randi(4);
            [x, y, v] = findNewLocation(divideLocation(1), divideLocation(2), plate, divideDir);
            if (v==1)
                nCells = nCells+1;
                plate(x,y) = 1;
                cellLocations(nCells,:) = [x,y];
            end
        end
    end
    set(im,'CData', plate);
    pause(.1);
end

end

function [x,y, valid] = findNewLocation(xin, yin, plate, direction)   
    x = xin;
    y = yin;
    valid = 1;
    % keep looking for new spot if current spot is occupied
    while( plate(x, y) == 1)
       switch direction
            case 1 % divide up
                y = y-1;
            case 2 % divide down
                y = y+1;
            case 3 % divide left
                x = x-1;
            case 4 % divide down
                x = x+1;
            otherwise
            warning('Invalid direction')
            x = xin;
            y = yin;
        return;
       end

       %if there has been a collision with a wall then just quit
       if y==0 || y==size(plate,2)+1 || x==0 || x==size(plate,1)+1 % hit the top
           x = xin; %return original values to say no division happend
           y = yin;
           valid = 0;
           return;
       end

    end


end

Note: Instead of thinking of pushing cells, I coded this in a way that leaves cells where they currently are and creates the new cell at the end of the row/column. Semantically its different but logically it has the same end result, as long as you don't care about the generations.

Answer 2

Inspired by an another question, I though of using image processing techniques to implement this simulation. Specifically we can use morphological dilation to spread the cancerous cells.

The idea is to dilate each pixel using a structuring element that looks like:

1 0 0
0 1 0
0 0 0

where the center is fixed, and the other 1 is placed at random at one of the other eight remaining positions. This would effectively extend the pixel in that direction.

The way the dilation is performed is by created a blank image, with only one pixel set, then accumulating all the results using a simple OR operation.

To speed things up, we don't need to consider every pixel, only those on the perimeter of the current blocks formed by the clusters of cancerous cells. The pixels on the inside are already surrounded by cancer cells, and would have no effect if dilated.

To speed even further, we perform the dilation on all pixels that are chosen to be extended in the same direction in one call. Thus every iteration, we perform at most 8 dilation operations.

This made the code relatively fast (I tested up to 1000x1000 grid). Also it maintains the same timing across all iterations (will not slow down as the grid starts to fill up).

Here is my implementation:

%# initial grid
img = false(500,500);

%# pick 10 random cells, and set them as cancerous
img(randi(numel(img),[10 1])) = true;

%# show initial image
hImg = imshow(img, 'Border','tight', 'InitialMag',100);

%# build all possible structing elements
%# each one dilates in one of the 8 possible directions
SE = repmat([0 0 0; 0 1 0; 0 0 0],[1 1 8]);
SE([1:4 6:9] + 9*(0:7)) = 1;

%# run simulation until all cells have cancer
BW = false(size(img));
while ~all(img(:)) && ishandle(hImg)
    %# find pixels on the perimeter of all "blocks"
    on = find(bwperim(img,8));

    %# percentage chance of division
    on = on( rand(size(on)) > 0.5 );    %# 50% probability of cell division
    if isempty(on), continue; end

    %# decide on a direction for each pixel
    d = randi(size(SE,3),[numel(on) 1]);

    %# group pixels according to direction chosen
    dd = accumarray(d, on, [8 1], @(x){x});

    %# dilate each group of pixels in the chosen directions
    %# to speed up, we perform one dilation for all pixels with same direction
    for i=1:8
        %# start with an image with only those pixels set
        BW(:) = false;
        BW(dd{i}) = true;

        %# dilate in the specified direction
        BW = imdilate(BW, SE(:,:,i));

        %# add results to final image
        img = img | BW;
    end

    %# show new image
    set(hImg, 'CData',img)
    drawnow
end

I also created an animation of the simulation on a 500x500 grid, with 10 random initial cancer cells (warning: the .gif image is approximately 1MB in size, so may take some time to load depending on your connection)