%************************************************************************** %***********UNIVERSITY OF ABERDEEN - iGEM CHEMOTAXIS MODEL 2009************ % %This is a simplified chemotaxis model for E. Coli in two dimensions. %This script can either create a single static image as a result, or can %produce an animated GIF, depending on the users needs. % %The data used has come from a variety of sources, listed in the references %section. The function describing how the tumbling frequency varies as the %bacteria move up a chemotactic gradient has been chosen without direct %scientific reference, other than the fact that the frequency varies %between 1 hertz and almost no tumbles per unit time. Although after seeing %bacteria chemotax under the microscope, many behave eratically and tumble %despite moving positively in the concentration gradient. % %The average speed of 27 micrometres per second is attained by randomly %varying the speed between 25 and 29 micrometres per second. This is based %on literature I have seen in which the speed does vary, and is not the %same for all bacteria, even in a relatively uniform culture. % %This program, despite not being agent based can account for collisions. %The collisions however are not in any way representative of real world %collisions, as there is no energy transferal or impact, simply if one %bacteria occupies another's space it will move away in a random direction %so that bacteria may not congregate on such a small point that the %concentration exceeds physical possibility. The bacteria are modelled by %spheres of a diameter of 2 micrometres of constant size. We decided not to %make an in depth collision model, as the main aim of the chemotaxis model %is to examine the behaviour of groups as they congregate around a source %of chemoattractant, and not of a tight biofilm. In order to examine how %useful the collision detection was, I recorded how many times collisions %occurred, and it was very very infrequently. Running the collision script %will take a lot longer, and the results are almost exactly the same. % %Exact data on the diffusion coefficient of Aspartate was not easy to come %by, however according to source [5] the diffusion constants of the amino %acids are very similar and fit a correlation to molecular mass. The %diffusion coefficient was estimated as 8*10^-10 m^2s^-1 from the graph. We %have not accounted for eddy diffusion, but it wouldn't be incorrect to %assume in many cases that diffusion occurs many hundreds of times faster %than this due to the bacteria mixing the fluid, and due to convection %currents, and agitation in the liquid. The value of c(0) in the diffusion %equation is the saturation concentration of aspartate in terms of moles %per metre cubed. % %Three iterations are equal to one second, and the bacteria check the %concentration every second. % %The data is stored in matrices of x and y coordinates where in pathx1(x,t) %x refers to the bacteria in question and t refers to the time. % %Please feel free to alter the program and play with it, the user input %mechanism has simply been put in place so that it may be used by people %with little or no knowledge of MATLAB. % %If making multiple gifs, you would be best to alter the name of the output %gif file in the imwrite command. % %In the static image output each green lightening like path, is the path %that the bacteria takes. % %REFERENCES: % %[1] Julius Adler. “Chemotaxis in Bacteria”. Annu. Rev. Biochem. 1975.44:341-356. % %[2] Nikhil Mittal, Elena O. Budrene, Michael P. Brenner, and Alexander van Oudenaarden. %“Motility of Escherichia coli cells in clusters formed by chemotactic aggregation”. %PNAS. Vol 100. no 23. Nov 11, 2003. % %[3] Kayo Maeda, Yasuo Imae, Jun-Ichi Shioi, and Fumio Oosawa. %“Effect of Temperature on Motility and Chemotaxis of Escherichia coli”. %JOURNAL OF BACTERIOLOGY, Vol. 127, No. 3 Sept. 1976, p. 1039-1046 % %[4] Jeffrey E. Segall, Steven M. Block, and Howard C. Berg. %“Temporal comparisons in bacterial chemotaxis”. Proc. Nati. Acad. %Sci. USA Biophysics Vol. 83, pp. 8987-8991, December 1986 % %[5] Alfred Polson. “On The Diffusion Constants of the Amino-Acids”. %Institute of Physical Chemistry, The University of Upsala, Sweden. 20 August 1937. % %**************Coded, Hacked, Written by Rory MacDonald******************** %************************************************************************** clear memory clear all %here are the user inputs: time = input('time in seconds that bacteria move ='); n = input('number of bacteria ='); collision = input('Do you want to inlcude collisions? (it will take longer to run) 1 = yes, 0 = no'); video = input('Do you want an animated GIF output? 1 = yes, 0 = no'); circle = input('Do you want to include a plot of where the signal chemical is enough for chemotaxis? 1 = yes, 0 = no'); speed = input('Do you want to have realistic slow speed, or increased speed? 1 = realistic, 0 = increased'); smile = input('Do you want the bacteria to be placed in an evenly distributed fashion, or in a smiley face? 1 = smile, 0 = even'); time2 = input('How long in seconds do you want to have had aspartate diffusing for before the simulation?'); %here is the declaration of variables and arrays steps = time*3 ; %3 fps - as in each iteration, or step, is a third of a second if speed == 1 div = 3; %dividing the speed in seconds, by the number of iterations per second puts the speed as it should be else div = 1; %dividing only by one gives a speed three times greater than in literature end pathx1 = ones(n,steps); pathy1 = ones(n,steps); dist1 = ones(n,steps); p1 = ones(n,steps); conc1 = ones(n,steps); D = 8.0*(10^-10); for x = 1:n %this piece of code assigns where each bacteria goes initially if smile == 1 cool(n) = rand; theta1 = -3 + 6*rand; if rand > 0.25 pathx1(x,1) = -8000 + cool(n)*16000; pathy1(x,1) = 4900*cos(cool(n)*pi - 5000); theta1 = -3 + 6*rand; else if rand > 0.5 pathx1(x,1) = -5000 + rand*100; pathy1(x,1) = 5000 + rand*100; theta1 = -3 + 6*rand; else pathx1(x,1) = 5000 + rand*100; pathy1(x,1) = 5000 + rand*100; theta1 = -3 + 6*rand; end end else pathx1(x,1) = -7500 + rand*15000; pathy1(x,1) = -7500 + rand*15000; theta1 = -3 + 6*rand; end for t = 2:steps %speed v = (25 + rand*4)/div; %changing the coordinates depending on direction and speed pathx1(x,t) = pathx1(x,t-1) + v*sin(theta1); pathy1(x,t) = pathy1(x,t-1) + v*cos(theta1); %measuring the distance from the source, in this case the origin dist1(x,t) = sqrt(((pathx1(x,t))^2)+((pathy1(x,t)^2))); %calculating the concentration at the source conc1(x,t) = (6938*(erfc((dist1(x,t)*10^-6)/(2*sqrt(D*((t+time2)/3)))))); %this part is the collision detection code. (it has to run once for all %bacteria before the current, and a second time for all after, so that %it does not detect a collision with itself if collision == 1 for z = 1:x-1 if pathx1(x,t) + 1 > pathx1(z,t) && pathy1(x,t) + 1 > pathy1(z,t) && pathx1(x,t) - 1 < pathx1(z,t) && pathy1(x,t) - 1 < pathy1(z,t) theta1 = 2*pi*rand; v = (25 + rand*4)/div; end end for z = x+1:n if pathx1(x,t) + 1 > pathx1(z,t) && pathy1(x,t) + 1 > pathy1(z,t) && pathx1(x,t) - 1 < pathx1(z,t) && pathy1(x,t) - 1 < pathy1(z,t) theta1 = 2*pi*rand; v = (25 + rand*4)/div; end end end %this part here describes the motion as the bacteria move up or down a %chemotactic gradient p1(x,t) = conc1(x,t)/conc1(x,t-1); %ratio of current to previous concentration if mod(t,3) == 0 %this is so that this only runs every third iteration if conc1(x,t) < 10^-8 %this makes the bacteria tumble if the concentration isn't high enough to be felt theta1 = 2*pi*rand; end if p1(x,t)<1 theta1 = 2*pi*rand; elseif 0.6*p1(x,t) 10^-8 radius(t) = 8000 - r; break else radius(t) = 0; end end end circlex = zeros(500,steps); circley = zeros(500,steps); for t = 2:steps for z = 1:505 circlex(z,t) = radius(t)*sin(z*((2*pi)/500)); circley(z,t) = radius(t)*cos(z*((2*pi)/500)); end end if video == 1 figure axis([-8000 8000 -8000 8000]) plot(pathx1(:,2),pathy1(:,2),'g.'); hold on plot(circlex(:,2),circley(:,2),'k'); hold off f = getframe; [im,map] = rgb2ind(f.cdata,256); im(1,1,1,20) = 0; axis equal for k = 2:steps box on plot(pathx1(:,k),pathy1(:,k),'g.'); hold on plot(circlex(:,k),circley(:,k),'k'); hold off axis([-8000 8000 -8000 8000]) f = getframe; im(:,:,:,k) = rgb2ind(f.cdata,map); end imwrite(im,map,'AberdeeniGEM.gif','DelayTime',0.00000000000000000001,'LoopCount',inf) else figure for x = 1:n plot(pathx1(x,:),pathy1(x,:),'g'); hold on plot(circlex(:,:),circley(:,:),'k'); end end else if video == 1 figure axis([-8000 8000 -8000 8000]) plot(pathx1(:,2),pathy1(:,2),'g.'); hold off f = getframe; [im,map] = rgb2ind(f.cdata,256); im(1,1,1,20) = 0; axis equal for k = 2:steps box on plot(pathx1(:,k),pathy1(:,k),'g.'); hold off axis([-8000 8000 -8000 8000]) f = getframe; im(:,:,:,k) = rgb2ind(f.cdata,map); end imwrite(im,map,'AberdeeniGEM.gif','DelayTime',0.00000000000000000001,'LoopCount',inf) else figure for x = 1:n plot(pathx1(x,:),pathy1(x,:),'g'); hold on end end end