## 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

Initial setup and declaration of variables...

%% AeroTaxis Implementation
% 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



## MAE3407 Finite Element Lab

This is the first of five computer sessions for the course MAE3407 Aircraft Structures II delivered at Monash University. These sessions are intended to merge the theory covered in class with the practical implementation of the numerical algorithms underlying the Finite Element Method (FEM). Each computer session will therefore contain exercises dealing with both analytical and numerical issues of finite elements. You are expected to solve and implement these exercises in MATLAB, document your solutions, and present them in a written report at the end of Week 10. By working through the lab exercises you will gain some insight into the synthesis of physics, mathematics, and computer science that form the basis of modern simulation techniques.