%************************************************************************** %***********UNIVERSITY OF ABERDEEN - iGEM CHEMOTAXIS MODEL 2009************ % %This is a simplified chemotaxis model for E. Coli in three dimensions. %This script can either create a single static image as a result, or can %produce an animated GIF, depending on the users needs. In this case the %bacteria can be placed in three different ways, and in each case they %chemotax in a pipe towards a crack from which diffuses aspartate. % %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, unlike the 2d version, does not account for collisions as it %is too computationally intensive. % %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. % %note: if you would like to get more info on how the script works, the 2d %version is more extensively commented % %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 time = input('time in seconds that bacteria move ='); n = input('number of bacteria ='); video = input('Do you want an animated GIF output? 1 = yes, 0 = no'); speed = input('Do you want to have realistic slow speed, or increased speed? 1 = realistic, 0 = increased'); view1 = input('Do you want the view to be one 3d view, or two 2d views from orthogonal view points? 1 = 2 views, 0 = 1 view'); time2 = input('How long in seconds do you want to have had aspartate diffusing for before the simulation?'); loca = input('How would you like the bacteria placed? 2 = left, 1 = right, 0 = even'); 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 D = 8.0*(10^-10); steps = time*3; %angle (initial) pathx1 = ones(n,steps); pathy1 = ones(n,steps); pathz1 = ones(n,steps); p1 = ones(n,steps); xdist1 = ones(n,steps); pipe1 = ones(n,steps); dist1 = ones(n,steps); completion = 0; for x = 1:n %this piece of code simply places the bacteria in the pipe in the %specified location. theta1 = 2*pi*rand; zeta1 = 2*pi*rand; if loca == 2 pathx1(x,1) = -5000 + rand*100; b = rand; pathy1(x,1) = 900 + 900*rand*sin(2*pi*b); pathz1(x,1) = 900*rand*cos(2*pi*b); end if loca == 1 pathx1(x,1) = 4900 + 100*rand; b = rand; pathy1(x,1) = 900 + 900*rand*sin(2*pi*b); pathz1(x,1) = 900*rand*cos(2*pi*b); end if loca == 0 pathx1(x,1) = -6000 + 12000*rand; b = rand; pathy1(x,1) = 900 + 900*rand*sin(2*pi*b); pathz1(x,1) = 900*rand*cos(2*pi*b); end for t = 2:steps %speed v = (25 + rand*4)/div; %this piece here updates the position of the bacteria at each time step pathx1(x,t) = pathx1(x,t-1) + v*sin(theta1); pathy1(x,t) = pathy1(x,t-1) + v*cos(theta1); pathz1(x,t) = pathz1(x,t-1) + v*cos(zeta1); xdist1(x,t) = pathx1(x,t); %this bit of the loop defines the crack if xdist1(x,t) >= -300 && xdist1(x,t) <= 300 xdist1(x,t) = 0; end %this defines the pipe as a cylinder dist1(x,t) = sqrt(((xdist1(x,t))^2)+((pathy1(x,t)^2))+((pathz1(x,t)^2))); pipe1(x,t) = sqrt((((pathy1(x,t)-900)^2))+((pathz1(x,t)^2))); if pipe1(x,t)>900 pathy1(x,t) = 900 + 900*((pathy1(x,t-1)-900)/pipe1(x,t-1)); pathz1(x,t) = 900*(pathz1(x,t-1)/pipe1(x,t-1)); end %this checks to see if the chemoattractant is in high enough concentration to be felt conc(x,t) = 6938*(erfc((dist1(x,t)*10^-6)/(2*sqrt(D*(t+time2))))); %probability of the bacteria deciding to tumble p1(x,t) = conc(x,t)/conc(x,t-1); q1 = rand; if mod(t,3) == 0 if conc(x,t)<(10^-8) theta1 = 2*pi*rand; zeta1 = 2*pi*rand; end if p1(x,t)<1 theta1 = 2*pi*rand; zeta1 = 2*pi*rand; elseif (0.6*p1(x,t)+((p1(x,t)^4)/8))