8

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!

Amro
  • 123,847
  • 25
  • 243
  • 454
DRG
  • 241
  • 1
  • 6
  • 1
    If you have a cancer blob which is already 3x3 and the cell at the centre of the blob divides, what is the rule (or what are the rules) for placing the daughter cell ? – High Performance Mark Jun 22 '12 at 13:47
  • Sounds like a problem that would lend itself well to a recursive solution instead of loops. – slayton Jun 22 '12 at 13:59
  • @HighPerformanceMark your question is answered in the original question: _"... 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"_. – Eitan T Jun 22 '12 at 14:00
  • Good question... the cell at centre (call this position 5, counting from top!) stays where it is, and pushes a daughter cell to anywhere in cells 1 - 4 or 6 - 9. As there is already a cell there, that cell gets pushed outside of the 3 x 3 grid - so if daughter cell went to position 4, the cell is 4 might would move to a new spot outside the 3 x 3 grid – DRG Jun 22 '12 at 14:01
  • 2
    What happens if the blob is touching the edge of a grid, say the left hand side) and a cell in the center divides left? Does that cell get lost or does force a division to the right? – slayton Jun 22 '12 at 14:02
  • Recursive might be the ideal system - any idea how to code that up ? I'm afraid I'm not a great coder and just about can handle multiple loops etc – DRG Jun 22 '12 at 14:02
  • 1
    @slayton regarding recursion: while it may be easier to code than loops, it might be slow like hell for large grids. – Eitan T Jun 22 '12 at 14:03
  • The edge of the grid ? Hmm.. initially we can allow it get lost if outside the 5000 x 5000 matrix, though i'd be interested in solutions in 3d space later on... right now, losses outside grid are ok – DRG Jun 22 '12 at 14:04
  • Does every cancer cell divide at every time step ? – High Performance Mark Jun 22 '12 at 14:11
  • Yes, at the moment they all do - later on I'll add in a percentage chance of division but for now each step has 2^n cancer cells.... – DRG Jun 22 '12 at 14:18
  • 1
    Do you need to keep track of which generation each cell was spawned in? Or are they all equivalent? If they are all equivalent, you can just place new daughter cells at the edge of the blob directly. Otherwise, it does get a bit more complicated. – tmpearce Jun 22 '12 at 14:30
  • +1 @tmpearce: that's exactly what I was about to ask. – Eitan T Jun 22 '12 at 14:31
  • At the moment generation isn't important. However, in later models it will be, as cells have different cycle times. Ugly, I know. Sorry. But I'd be happy with any outline method!! – DRG Jun 22 '12 at 14:35
  • When a daugher cell pushes an existing cancer cell, does it get pushed to a random adjacent cell? – Eitan T Jun 22 '12 at 14:43
  • Yes indeed - the "bumping" cell takes place of occupied cell in a domino effect – DRG Jun 22 '12 at 14:47

2 Answers2

4

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.

slayton
  • 20,123
  • 10
  • 60
  • 89
  • Seems fantastic - just trying to get it to run now - pardon my ignorance, but do I have to define plateSize as the grid (ie: 5000) and what is pDivide ? – DRG Jun 22 '12 at 14:48
  • If you want a grid of 5000x5000 then `plateSize` is 5000. `pDivide` is the probability that a given cell will divide each generation. If you set it to 1 then each cell will always divide. If you set it to .1 then each cell will have a 10% chance of dividing each generation. – slayton Jun 22 '12 at 14:52
  • Cheers - it's going really well, except it keeps telling me there's not enough input arguments on line " if divide >= pDivide" any ideas ? – DRG Jun 22 '12 at 15:03
  • you need to pass it two input arguments ie. `simulateCancer(5000,1);` – slayton Jun 22 '12 at 15:18
  • I'm sorry, I really don't know what this means in practice - I essentially copy and pasted your code, added extra lines for pDivide = 1 and just for testing plateSize = 500; apologies for not getting it quicker, it keeps returning an input argument problem – DRG Jun 22 '12 at 15:24
  • 1
    create a file called simulateCancer.m, paste my code into it, and then execute the code with `simulateCancer(5000,1)` from matlab's command window – slayton Jun 22 '12 at 15:38
  • let us [continue this discussion in chat](http://chat.stackoverflow.com/rooms/12916/discussion-between-user1474967-and-slayton) – DRG Jun 22 '12 at 15:44
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)

simulation_animation

Community
  • 1
  • 1
Amro
  • 123,847
  • 25
  • 243
  • 454