Today, I’ll show how to simply visualize 3D tumor using MATLAB’s isosurface function combined with a simple blurring technique.

First, we will modify the code from the previous post in order to simulate tumor in 3D. There are only a few changes: redefinition of the domain and the neighborhood; generating permutations on the fly instead of storing them in *Pms* varaible; and the way in which the initial cell is placed.

N = 200; %cubic domain dimension nSteps = 30*24; %number of simulation steps pprol = 1/24; %probability of proliferating pmig = 10/24; %probability of migrating pdeath = 1/100; %probability of dying pmax = 10; %proliferation capacity ps = 3/10; %probability of symmetric division L = false(N,N,N); %domain definition L(1,:,:) = true; L(end,:,:) = true; %filling boundary L(:,1,:) = true; L(:,end,:) = true; %filling boundary L(:,:,1) = true; L(:,:,end) = true; %filling boundary L(N*N*round(N/2)+N*round(N/2)+round(N/2)) = true; cells = int32(N*N*round(N/2)+N*round(N/2)+round(N/2)); cellsIsStem = true; cellsPmax = uint8(pmax); aux = int32([[-N-1 -N -N+1 -1 1 N-1 N N+1] ... [-N-1 -N -N+1 0 -1 1 N-1 N N+1]-N*N ... [-N-1 -N -N+1 0 -1 1 N-1 N N+1]+N*N])'; %indices to heighborhood for i = 1:nSteps sh = randperm(length(cells)); cells = cells(sh); cellsIsStem = cellsIsStem(sh); cellsPmax = cellsPmax(sh); Pms = cell2mat(arrayfun(@(x)randperm(26),(1:length(cells))','UniformOutput',0))'; S = bsxfun(@plus,cells,aux(Pms)); S(L(S)) = 0; %setting occupied spots to 0 indxF = find(any(S)); %selecting cells with at least one free spot nC = length(indxF); %number of cells with free spot P = rand(1,nC)<pprol; %proliferation Ps = P & rand(1,nC)<ps & cellsIsStem(indxF);%symmetric division De = P & (cellsPmax(indxF) == 0);%proliferation exhaution D = P & (rand(1,nC)<pdeath) & ~cellsIsStem(indxF); %death at proliferation attempt M = ~P & (rand(1,nC)<pmig); %go when no grow del = D | De; %cells to delete act = find((P | M) & ~del); %indices to the cells that will do something for ii = act %only for those that will do anything ngh = S(:,indxF(ii)); ngh(ngh==0) = []; indO = find(~L(ngh),1,'first'); %selecting free spot if ~isempty(indO) %if there is still a free spot L(ngh(indO)) = true; if P(ii) %proliferation cells = [cells uint32(ngh(indO))]; if Ps(ii) %symmetric division cellsIsStem = [cellsIsStem true]; cellsPmax = [cellsPmax cellsPmax(indxF(ii))]; else cellsIsStem = [cellsIsStem false]; cellsPmax = [cellsPmax cellsPmax(indxF(ii))-1]; if ~cellsIsStem(indxF(ii)) cellsPmax(indxF(ii)) = cellsPmax(indxF(ii))-1; end end else %migration L(cells(indxF(ii))) = false; cells(indxF(ii)) = uint32(ngh(indO)); end end end if ~isempty(del) %updating death L(cells(indxF(del))) = false; cells(indxF(del)) = []; cellsIsStem(indxF(del)) = []; cellsPmax(indxF(del)) = []; end end

Using the above code I simulated two tumors that have different growth characteristics, both with pmig = 5/24 and one with ps = 3/10 and the other one with ps = 5/100. Now we need to visualize them.

Having the whole tumor stored in the *L* variable we can use isosurface function straight away and plot the tumor as it is. The only thing that we need to do is to remove the cells from the *L* boundary.

function visualize( N, L) %clearing cells from the boundary L(1,:,:) = false; L(end,:,:) = false; L(:,1,:) = false; L(:,end,:) = false; L(:,:,1) = false; L(:,:,end) = false; %calculating isosurfaces and plotting p = patch(isosurface(1:N,1:N,1:N,L,0.25)); isonormals(1:N,1:N,1:N,L,p) set(p,'FaceColor','red','EdgeColor','none'); xlim([1 N]); ylim([1 N]); zlim([1 N]); view([90 0]); camlight lighting gouraud end

The results of the above function are plotted below.

As you can see the visualization doesn’t really allow to see the 3D structure of the tumors, as individual separated cells are plotted and there is no clear 3D surface. We will fix that by using a simple blurring technique. The basic idea is that we will replace each site value by the average value from its neighborhood (several times in a loop). That will allow to delete separated cells and make 3D surface smoother. We add an input variable *blLevels* to the function in order to define how may iterations of that procedure should be performed.

function visualize( N, L, cells, blLevels ) %clearing cells from the boundary L(1,:,:) = false; L(end,:,:) = false; L(:,1,:) = false; L(:,end,:) = false; L(:,:,1) = false; L(:,:,end) = false; if blLevels %if perform basic blur %auxilary variable with indices to the cell negiborhood aux = int32([[-N-1 -N -N+1 -1 1 N-1 N N+1] ... [-N-1 -N -N+1 0 -1 1 N-1 N N+1]-N*N ... [-N-1 -N -N+1 0 -1 1 N-1 N N+1]+N*N])'; %creating indices to cells and their beignorhoods S = [cells unique(reshape(bsxfun(@plus,cells,aux),1,[]))]; %creating indices to neigborhood of indices in S S2 = bsxfun(@plus,S,aux); %making sure that indices are still within the lattice S2(S2<1) = []; S2(S2>N*N*N) = []; %changing lattice from logical variable to float L = single(L); for i = 1:blLevels %for number of blurs L(S) = mean(L(S2)); %taking the average of neighborhood end end %calculating isosurfaces and plotting p = patch(isosurface(1:N,1:N,1:N,L,0.25)); isonormals(1:N,1:N,1:N,L,p) set(p,'FaceColor','red','EdgeColor','none'); xlim([1 N]); ylim([1 N]); zlim([1 N]); view([90 0]); camlight lighting gouraud end

If we apply *blLevels* = 1 we get the following, nicer plot of the tumor.

For *blLevels* = 2 we get the following plots.

Finally, when applying three iterations of blur (*blLevels* = 3) we get 3D plots nicely showing the tumor shape.

What is also great it takes only about 3 seconds to plot a tumor with more than quarter of a million of cells (left one).

Reblogged this on CancerEvo and commented:

Nice post by Jan Poleszcuk about 3D visualisation in Matlab. Now if only we could persuade him to do this with Python or other open source tools…

LikeLike

Hey can you please give me an example for the inputs to visualize the data

LikeLike

You have the code that generates the input to the visualization function in the first part of the code, but in general the input variable” L” is a 3D logical matrix that have true values only in places occupied by cells + on the boundary and “cells” is an accompanying vector with positions of cells using linear indexing. Of course you can construct “cells” from “L” and vice versa if you want to reduce number of input parameters.

LikeLike