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;

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.

Click to download Lab 1 notes

Click to download Lab 2 notes

Click to download Lab 3 notes

Click to download Lab 4 notes

Click to download Lab 5 notes

Click to download getElement.m

Click to download lab3.mat

Click to download lab4.mat

Click to download calculate_SXX.m

Click to download createHolyRectangle.m

Diablo 3 On Debian Squeeze

I got sick of putting up with Windows just to play Diablo 3, so I looked into running it on Debian with Wine. It works great, and surprisingly doesn't have lag spikes like it did under Windows 7.The process below is probably the same or even simpler for Ubuntu users.

Firstly, you will need the latest version of Wine with the "AcceptEx fix" patchset by Erich Hoover. As of Wine 1.5.6 the patches submitted by Erich were commited. If you have any experience with running games with Wine, you'll know that certain games run better with particular versions of Wine. On a distribution like Debian, which has a package management system, that means you'll need multiple, parallel installations of Wine for each game. This can be confusing and messy so it's best to either compile Wine yourself in separate directories or have PlayOnLinux do it for you.

PlayOnLinux can be installed on Debian by

apt-get install playonlinux

PlayOnLinux has a large database of scripts and workthroughs for a large number of Windows games and applications. This means, you simply ask PlayOnLinux to install a game and it will hopefully do everything for you. This includes downloading the correct version of Wine and other addons, etc, etc. I won't describe how to use PlayOnLinux because it has a friendly GUI that you can work your way around.

I found that PlayOnLinux started the Launcher and downloaded Diablo 3 without any problems. However, as soon as it was ready to play it did nothing when I clicked on "Play". So I manually started it from the command line interface, and found the following error being repeated indefinitely:

Direct3D9 is not available without OpenGL.

It turns out this problem is caused by the lack of 32bit opengl libraries (I'm running 64bit debian), more specifically nvidia-glx-ia32 package. I tried numerous times to install this via the package management, but it failed miserably during the post configuration step. I got errors like:

Package libgl1-nvidia-glx-ia32 is not configured yet

After a bit of tinkering, I remembered I had installed cuda drivers for GPU number crunching. Maybe there was a clash? Anyway, I removed all debian nvidia packages, disabled X server on start up and restarted Debian. Using the Nvidia drivers from, I installed the latest modules and libraries. Interestingly, the installer does ask you if you want to install 32bit compatibility libraries. Make sure you click Yes 🙂

After getting X server up, Diablo 3 worked fine. After killing Diablo I can safely say it runs beautifully- no crashes or glitches.

Good luck!

Yay, it works
Yay, it works

Debian Installation Tips

Recently I've switched back to my first and favourite Linux distribution Debian. The latest release is called Debian Squeeze or Debian 6.0. So far it has been great. I have found that after setting up Debian on my server and multiple workstations the setup process isn't exactly complete when the installer is finished. So I've put together some hints/tips that I do on every Debian installation.

1. Software Sources

We need to modify the software sources to include the non-free and contrib repositories. You can do this directly via CLI or using the Software Source app under System -> Administration. Make sure you check the two boxes "DFSG-compatible Software with Non-free Dependencies (contrib)" and "Non-DFSG-compatible Software (non-free)". If you have installed Debian using a DVD you should untick the cdrom options under the Third-Part-Software tab. Finally, update the database manually if it doesn't automatically:

apt-get update

2. Nvidia Drivers

The best reference for this would be the wiki available here. I've listed a cut down version of the required commands for a typical setup below for your copy and pasting pleasure:

apt-get install module-assistant nvidia-kernel-common

m-a auto-install nvidia-kernel${VERSION}-source

apt-get install nvidia-glx${VERSION} nvidia-xconfig

Then simply run nvidia-xconfig and restart X server.

3. Fonts

You will notice that Debian doesn't ship with the pretty fonts. This is easily noticable when browsing the Web. To install fonts:

apt-get install ttf-mscorefonts-installer

4. 32 bit 3D libraries

If you have installed Debian 64bit like I always do, then you will surely have a problem with either Flash or other applications such as Google Earth. For example, you might be able to use google maps but Google street view simply doesn't work. Or you might see that Google Earth shows only a blank screen. Since you are running 64bit Debian you need to install the 32bit compatibility libraries:

apt-get install nvidia-glx-ia32 ia32-libs-gtk lib32nss-mdns ia32-libs lib32ncurses5

How to calculate Optical Forces using MEEP

The optical force acting on an arbitrary particle in an arbitrary incident field can be calculated using the integral

 \mathbf{F} = \int_{S} \left\lbrace \mathbf{T}(\mathbf{r},t) \right\rbrace \cdot \mathbf{n}(\mathbf{r}) dS

where  \left\lbrace \cdot \right\rbrace ,  S and  \mathbf{T} represent the time-averaged operator, surface enclosing the particle and Maxwell's stress tensor, respectively. Maxwell's stress tensor is given as

 \mathbf{T} = \left[ \varepsilon_0 \mathbf{E}\mathbf{E} - \mu_0 \mathbf{H}\mathbf{H} - \frac{1}{2}(\varepsilon_0 E^{2} + \mu_0 H^2) \mathbf{I} \right]

which requires all field components to be calculated. You can obtain the solved field components in a number of ways. We can use the data from MEEP to calculate the stress tensor and hence the optical forces. If we are interested in finding the forces at a grid of points then we must run the MEEP simulation for each point.

I have provided a set of files which calculates the forces a particle experiences in a standing wave. It should be fairly easy to modify to suit your simulation.

- : this file runs everything and plots the forces. You will need numpy/scipy installed.

- standing.ctl : meep simulation file that generates stress tensor for a particle at different locations (you can see txx, txy, etc definitions)

- standing-full.ctl : meep simulation to obtain field density (only for plotting)

- h5tonumpy : convert h5 dataset to numpy file.

This method of obtaining optical forces is very slow and there are inherent inaccuracies. There was a paper that mentioned this. I will provide I link to it when I find it again. Essentially we are only using FDTD to obtain the fields which can be obtained using other scattering codes/methods. However, for more complicated geometries and/or materials (such as plasmonics/meta-materials) the FDTD brute force method can be useful.

Please note: The latest version of MEEP now has an inbuilt function to calculate Maxwell's stress tensor. It is much easier to calculate forces using that. However I leave this here for legacy reasons.

Download standing-wave-example-tar


I have been asked by numerous people if it is possible to implement a traffic quota and shaping system at home. It seems we all know of or have had direct experience with bandwidth hogs who use up all the monthly bandwidth downloading their favourite Anime, leaving the rest of us with shaped internet. This howto aims to provide a set of instructions and tools to deploy a standalone Linux based Traffic accounting and Shaping System.

Download ZOIDtraf reference PDF

Download ZOIDtraf source