%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% E. coli Bacterial Swarm Foraging for Optimization
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Kevin Passino
% Version: 6/1/00
%
% This program simulates the minimization of a simple function with
% chemotaxis, swarming, reproduction, and elimination/dispersal of a
% E. coli bacterial population.
%
% To change the nutrientsfunc search on it. For
% example, change it to nutrientsfunc1 to study another type of swarm behavior.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all % Initialize memory
p=2; % Dimension of the search space
S=50; % The number of bacteria in the population (for convenience, require this to be an
% an even number)
Nc=100; % Number of chemotactic steps per bacteria lifetime (between reproduction steps), assumed
% for convenience to be the same for every bacteria in the population
Ns=4; % Limits the length of a swim when it is on a gradient
Nre=4; % The number of reproduction steps (right now the plotting is designed for Nre=4)
Sr=S/2; % The number of bacteria reproductions (splits) per generation (this
% choice keeps the number of bacteria constant)
Ned=2; % The number of elimination-dispersal events (Nre reproduction steps in between each event)
ped=0.25; % The probabilty that each bacteria will be eliminated/dispersed (assume that
% elimination/dipersal events occur at a frequency such that there can be
% several generations of bacteria before an elimination/dispersal event but
% for convenience make the elimination/dispersal events occur immediately after
% reproduction)
flag=2; % If flag=0 indicates that will have nutrients and cell-cell attraction
% If flag=1 indicates that will have no (zero) nutrients and only cell-cell attraction
% If flag=2 indicates that will have nutrients and no cell-cell attraction
% Initial population
P(:,:,:,:,:)=0*ones(p,S,Nc,Nre,Ned); % First, allocate needed memory
% Initialize locations of bacteria all at the center (for studying one swarming case - when nutrientsfunc1 is used)
%for m=1:S
% P(:,m,1,1,1)=[15;15];
%end
% Another initialization possibility: Randomly place on domain:
for m=1:S
P(:,m,1,1,1)=(15*((2*round(rand(p,1))-1).*rand(p,1))+[15;15]);
end
% Next, initialize the parameters of the bacteria that govern
% part of the chemotactic behavior
C=0*ones(S,Nre); % Allocate memory
% Set the basic run length step (one step is taken of this size if it does not go up a gradient
% but if it does go up then it can take as many as Ns such steps). C(i,k) is the step size for
% the ith bacteria at the kth reproduction step. For now, the step size is assumed to be constant since
% we assume that perfect copies of bacteria are made, but later you can add evolutionary effects to modify
% this and perhaps Ns and Nc.
runlengthunit=0.1;
C(:,1)=runlengthunit*ones(S,1);
% Allocate memory for cost function:
J=0*ones(S,Nc,Nre,Ned);
Jhealth=0*ones(S,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%---------------------------------
% Elimination-dispersal loop:
%---------------------------------
for ell=1:Ned
%---------------------------------
% Reproduction loop:
%---------------------------------
for k=1:Nre
%---------------------------------
% Swim/tumble (chemotaxis) loop:
%---------------------------------
for j=1:Nc
for i=1:S % For each bacterium
% Compute the nutrient concentration at the current location of each bacterium
J(i,j,k,ell)=nutrientsfunc(P(:,i,j,k,ell),flag);
% J(i,j,k,ell)=nutrientsfunc1(P(:,i,j,k,ell),flag);
% Next, add on cell-cell attraction effect:
J(i,j,k,ell)=J(i,j,k,ell)+bact_cellcell_attract_func(P(:,i,j,k,ell),P(:,:,j,k,ell),S,flag);
%-----------
% Tumble:
%-----------
Jlast=J(i,j,k,ell); % Initialize the nutrient concentration to be the one at the tumble
% (to be used below when test if going up gradient so a run should take place)
% First, generate a random direction
Delta(:,i)=(2*round(rand(p,1))-1).*rand(p,1);
% Next, move all the bacteria by a small amount in the direction that the tumble resulted in
% (this implements the "searching" behavior in a homogeneous medium)
P(:,i,j+1,k,ell)=P(:,i,j,k,ell)+C(i,k)*Delta(:,i)/sqrt(Delta(:,i)'*Delta(:,i));
% This adds a unit vector in the random direction, scaled
% by the step size C(i,k)
%---------------------------------------------------------------------
% Swim (for bacteria that seem to be headed in the right direction):
%---------------------------------------------------------------------
J(i,j+1,k,ell)=nutrientsfunc(P(:,i,j+1,k,ell),flag); % Nutrient concentration for each bacterium after
% J(i,j+1,k,ell)=nutrientsfunc1(P(:,i,j+1,k,ell),flag); % Nutrient concentration for each bacterium after
% a small step (used by the bacterium to
% decide if it should keep swimming)
% Next, add on cell-cell attraction effect:
J(i,j+1,k,ell)=J(i,j+1,k,ell)+bact_cellcell_attract_func(P(:,i,j+1,k,ell),P(:,:,j+1,k,ell),S,flag);
m=0; % Initialize counter for swim length
while mrand % Generate random number and if ped bigger than it then eliminate/disperse
P(:,m,1,1,ell+1)=(15*((2*round(rand(p,1))-1).*rand(p,1))+[15;15]);
else
P(:,m,1,1,ell+1)=P(:,m,1,Nre+1,ell); % Bacteria that are not dispersed
end
end
%---------------------------------
end % ell=1:Ned
%---------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot the function we are seeking the minimum of:
x=0:31/100:30; % For our function the range of values we are considering
y=x;
% Compute the function that we are trying to find the minimum of.
for jj=1:length(x)
for ii=1:length(y)
z(ii,jj)=nutrientsfunc([x(jj);y(ii)],flag);
% z(ii,jj)=nutrientsfunc1([x(jj);y(ii)],flag);
end
end
% First, show the actual function to be maximized
figure(1)
clf
surf(x,y,z);
colormap(jet)
% Use next line for generating plots to put in black and white documents.
colormap(white);
xlabel('x=\theta_1');
ylabel('y=\theta_2');
zlabel('z=J');
title('Nutrient concentration (valleys=food, peaks=noxious)');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Next, provide some plots of the results of the simulation.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
t=0:Nc; % For use in plotting (makes t=0 correspond to the 1 index and plots to Nc+1)
% As is Figure 2 (4) shows parameter trajectories for 4 generations, then there is an elimination/dispersal event
% Figure 3 (5) shows parameter trajectories for the following 4 generations (after the elim/disp event)
for kk=1:Ned
figure(kk+1)
clf
for mm=1:Nre
subplot(2,2,mm)
for nn=1:S % Plot all bacteria trajectories for generation mm
plot(t,squeeze(P(1,nn,:,mm,kk)),t,squeeze(P(2,nn,:,mm,kk)))
%plot(t,squeeze(P(1,nn,:,mm,kk)),'-',t,squeeze(P(2,nn,:,mm,kk)),'-')
axis([min(t) max(t) 0 30])
hold on
end
T=num2str(mm);
T=strcat('Bacteria trajectories, Generation=',T);
title(T)
xlabel('Iteration, j')
ylabel('\theta_1, \theta_2')
hold off
end
end
for kk=1:Ned
figure(Ned+kk+1)
clf
for mm=1:Nre
subplot(2,2,mm)
contour(x,y,z,25)
colormap(jet)
for nn=1:S % Plot all bacteria trajectories for generation mm
hold on
plot(squeeze(P(1,nn,:,mm,kk)),squeeze(P(2,nn,:,mm,kk)))
%plot(squeeze(P(1,nn,:,mm,kk)),squeeze(P(2,nn,:,mm,kk)),'-')
axis([0 30 0 30])
end
T=num2str(mm);
T=strcat('Bacteria trajectories, Generation=',T);
title(T)
% Use next line for generating plots to put in black and white documents.
%colormap(gray);
xlabel('\theta_1');
ylabel('\theta_2');
hold off
end
end
%%%%%%%%%%%%
pause % Can leave this in if want to avoid movie (then hit control-C)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Next, show a movie of the chemotactic steps:
figure(2*Ned+2)
clf
contour(x,y,z,25)
colormap(jet)
axis([0,30,0,30]);
xlabel('\theta_1');
ylabel('\theta_2');
title('Bacteria movements');
hold on
M = moviein(Nc);
for j=1:Nc;
% Can change generation step and elimination-dispersal step on next line.
% Currently for 1,1 - the first generation in the first elimination dispersal step
for i=1:S
v=plot(squeeze(P(1,i,j:j+1,1,1)),squeeze(P(2,i,j:j+1,1,1)),'-');
set(v,'MarkerSize',3);
end
M(:,j)=getframe;
end;
%movie(M,0)
%save bacteria_swarm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% End of program
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%