Mechanics of bacterial aerotaxis and its simulation

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

 \frac{dl}{dt} = v \frac{dl}{dx} + f_{rl} r - f_{lr} l

 \frac{dr}{dt} = -v \frac{dr}{dx} + f_{lr} l - f_{rl} r

Eq. 1

where v=40 \mu m/s is the constant swimming speed of the bacteria, f_{rl} and f_{lr} 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 c_{min} = 0.015 , lower detectable oxygen concentration \hat{c}_{min} = 0.0005 , upper favorable oxygen concentration c_{max} = 0.025 , and upper favorable oxygen concentration \hat{c}_{max} = 0.25.

Consequently, the reversal energy can be described using the lower reversal frequency f=0.1s^{-1} and the upper reversal frequency F=0.5 s^{-1} via

where f(E,dE/dt) = F when dE/dt < 0 and f(E,dE/dt) when dE/dt>0. In the process, the bacterial cells contribute to the balance of oxygen concentration L in the capillary via

\frac{dL}{dt} = D \frac{d^2 L}{dx^2} - \kappa u(L) b

Eq. 2

where D=2\times 10^{3} \mu m/s^2 is the diffusion coefficient of oxygen in water, and \kappa = 10^{-9}(\mu M ml)/(s \dot cell) the oxygen consumption rate of the cells. The step function introduced operates such that u(L) = 1 at c>0, and u(L) = 0 for c\leq 0.

Let's discretize the spatial domain x and time domain t into n_x and n_y points with a step size \Delta x and \Delta t, 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,

 r_{j}^{n+1} = r_{j}^{n} - \frac{1}{2} K_{r} \left( r_{j+1}^{n} - r_{j-1}^{n} \right) + \frac{1}{2} K_{r}^{2} \left( r_{j+1}^{n} - 2r_{j}^{n} + r_{j-1}^{n} \right) - \Delta t f_{rl} r_{j}^{n} + \Delta t f_{lr} l_{j}^{n}

 l_{j}^{n+1} = l_{j}^{n} - \frac{1}{2} K_{l} \left( l_{j+1}^{n} - l_{j-1}^{n} \right) + \frac{1}{2} K_{l}^{2} \left( l_{j+1}^{n} - 2l_{j}^{n} + l_{j-1}^{n} \right) - \Delta t f_{lr} l_{j}^{n} + \Delta t f_{rl} r_{j}^{n}

Eq. 3

where the superscripts and subscripts represent the nth time step and the jth spatial step, respectively. The parameter K_{r} = -K_{l} = -v \Delta t / \Delta x can be used to define the Courant-Friedrichs-Level (CFL) condition for stability, which requires that \| K \| \leq 1. We find the aspect ratio of  \Delta t = \Delta x / v 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

 L_{j}^{n} = L_{j}^{n+1} \left( 1 + 2 \sigma \right) - \sigma \left( L_{j+1}^{n+1} + L_{j-1}^{n+1} \right) + \kappa \Delta t u(L_{j}^{n+1}) \left( r_{j}^{n+1} + l_{j}^{n+1} \right)

Eq. 4

where  \sigma = D \Delta t / (\Delta x)^2. The Neumann boundary condition  \frac{dL}{dx} = \gamma on the right end of the capillary tube can be described using a ghost-point n_x + 1 which gives L_{n_x+1}^{n+1}= 2 \Delta x \gamma + L_{n_x-1}^{n+1}. Substituting this into Eq. (4) when evaluated such that j=n_x and eliminating L_{n_x + 1}^{n+1}, reveals two additional components -2\Delta x \gamma - \sigma L_{n_x-1}^{n+1} 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 \gamma = 0.

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,


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

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
    L(end,n) = R(end,n);
    R(end,n) = 0.0;
    R(1,n) = L(1,n);
    L(1,n) = 0.0;



and finally, plot and view

close all
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])
    i = i + 1;