In an assay where bacteria is attracted to a specific low concentration of oxygen, the gradient of oxygen that is established will allow aerotaxis to guide the bacteria towards a concentration of oxygen that is optimal for energy production. This behaviour is governed by the thesis that the bacteria that leaves the aerotactic band through the high-oxygen or low-oxygen boundary will experience a loss of energy, which then increases the reversal frequency resulting in a return to this band. Upon reversal into the band, the energy increase causes a suppression of this reversal. Consequently, this leads to the trapping of bacteria within the band.

In a 1D model, we assume gas infusion is from the left to the right end of the capillary. Hence, the respective bacterial densities *l*(*x*, *t*) and *r*(*x*, *t*) can be given by

Eq. 1

where is the constant swimming speed of the bacteria, and are the frequencies of the reversals from right to left and from left to right. These reversals depend on the local oxygen concentration and its derivatives. It can be assumed that the energy E of the bacteria cells range between Emin and Emax, according to values of the lower favorable oxygen concentration , lower detectable oxygen concentration , upper favorable oxygen concentration , and upper favorable oxygen concentration .

Consequently, the reversal energy can be described using the lower reversal frequency and the upper reversal frequency via

where when and when . In the process, the bacterial cells contribute to the balance of oxygen concentration in the capillary via

Eq. 2

where is the diffusion coefficient of oxygen in water, and the oxygen consumption rate of the cells. The step function introduced operates such that at , and for .

Let's discretize the spatial domain and time domain into and points with a step size and , respectively. Using the Lax-Wendroff finite difference scheme, which is second-order accurate in both space and time, we obtain from Eq. (1) an explicit numerical scheme,

Eq. 3

where the superscripts and subscripts represent the nth time step and the jth spatial step, respectively. The parameter can be used to define the Courant-Friedrichs-Level (CFL) condition for stability, which requires that . We find the aspect ratio of suitably stable and fast. To conserve the total bacteria population we apply a boundary condition at the end points such that the right going population switches to left going at the right-most boundary, and vice versa. An implicit second-order accurate numerical scheme is applied to Eq. (2) obtaining

Eq. 4

where . The Neumann boundary condition on the right end of the capillary tube can be described using a ghost-point which gives . Substituting this into Eq. (4) when evaluated such that and eliminating , reveals two additional components that need to be added to the right hand side of Eq. (4) to satisfy the Neumann boundary condition. To physically reflect the no-flux condition at the capillary tube end, we set .

**Matlab Implementation**

Download AeroTaxis MATLAB Implementation

If you find this code useful, please cite this article..

Initial setup and declaration of variables...

%% AeroTaxis Implementation % by Murat Muradoglu, murat.neyzen@gmail.com % www.torrkish.com % http://www.torrkish.com/index.php/2016/03/29/mechanics-of-bacterial-aerotaxis-and-its-simulation/ clear delX = 1/200; % Bacteria swimming speed v = 1.0; % Diffusion Constant D = 0.5; % Courant-Friedrichs-Levy (CFL) condition % states that C = |v|*(delT/delX) <= 1 % but we find that it is only reliably stable % when C = 1 . delT = delX./v; K1 = -v*delT/(delX); K2 = v*delT/(delX); sigma = (D*delT)/(delX.^2); F = 1.25; % UPPER reversal frequency f = 0.25; % LOWER reversal frequency Kap = 1.0; % Oxygen consumption rate by bateria t = 0:delT:15; x= 0:delX:10; tsteps = length(t); xsteps = length(x); R = zeros(xsteps,tsteps); % right moving density L = zeros(xsteps,tsteps); % left moving OX = zeros(xsteps,tsteps); % oxygen bE2 = zeros(xsteps,tsteps); % Energy calcs. % Oxygen initial conditions % Dirichlet boundary condition OX(:,1) = 0; OX(end,1) = 1; % Right and Left population initial conditions R(:,1) = 0.1; L(:,1) = 0.1; % Construct matrix Oxygen-diffusion A1 = diag(ones(xsteps-1,1),-1)*-1*sigma; A2 = diag(ones(xsteps,1))*(1+2*sigma); A3 = diag(ones(xsteps-1,1),1)*-1*sigma; A = sparse(A1 + A2 + A3); % Boundary Conditions on Oxygen-diffusion matrix A(end,:) = 0; A(end,end) = 1; % Oxygen-Neumann condition A(1,2) = A(1,2)*2; % Store inverse of A invA = inv(A); % Eating matrix E = sparse(diag(ones(xsteps-1,1),1)*delT*Kap); % right propagator B1 = diag(ones(xsteps-1,1),1)*0.5*K1*(1+K1); B2 = diag(ones(xsteps,1))*(1-K1.^2); B3 = diag(ones(xsteps-1,1),-1)*-0.5*K1*(1-K1); B = sparse(B1 + B2 + B3); % left propagator C1 = diag(ones(xsteps-1,1),1)*0.5*K2*(1+K2); C2 = diag(ones(xsteps,1))*(1-K2.^2); C3 = diag(ones(xsteps-1,1),-1)*-0.5*K2*(1-K2); C = sparse(C1 + C2 + C3);

The loop. Runs rather quick. ~3 seconds on my laptop

tic for n=2:tsteps % Propagate Bacteria R(:,n) = B*R(:,n-1); L(:,n) = C*L(:,n-1); % Calculate energy experienced by bacteria bE2(:,n) = calcEnergy(OX(:,n-1),x); % Calculate respective reversal functions [frl,flr] = calcReversal(x,bE2(:,n),F,f); % Loss from right->left and left-> right Rloss = frl.*R(:,n); Lloss = flr.*L(:,n); % Apply reversals.. R(:,n) = R(:,n) - delT.*Rloss + delT.*Lloss; L(:,n) = L(:,n) - delT.*Lloss + delT.*Rloss; % Total bacteria population Total = R(:,n) + L(:,n); % Diffuse Oxygen % forward eating OX(:,n) = A\OX(:,n-1); % Much faster to do A\Ox than invA*Ox, even when invA is stored inverse OmgTotal = (OX(:,n) > 0).*Total; % Eat only if there is enough oxygen OX(:,n) = OX(:,n) - E*OmgTotal; OX(OX(:,n) < 0,n) = 0; % Can't eat "negative" oxygen % Boundary Condition %R(end,n) = R(1,n); % periodicity Boundary condition % RIGHT WALL L(end,n) = R(end,n); R(end,n) = 0.0; % LEFT WALL R(1,n) = L(1,n); L(1,n) = 0.0; end toc

and finally, plot and view

close all figure i = 1; for p=1:10:tsteps plot(x, [R(:,p) + L(:,p) OX(:,p) ]) Totals(i) = sum(R(:,p) + L(:,p)); %% Total should converge to a certain value %fprintf('Total: %d\n',Totals(i)); xlim([0 10]) drawnow pause(0.1) i = i + 1; end