A simple model for chemotaxis

Learning Scientific Programming with Python (2nd edition)

P6.6.3: A simple model for chemotaxis

Question P6.6.3

Some bacteria, such as E. coli, possess helical flagella which enable them to move towards attractants such as nutrients, process known as chemotaxis. When the flagella rotate anti-clockwise the bacterium is propelled forwards; when they rotate clockwise, it tumbles randomly, changing its orientation. A combination of such movements enables the bacterium to perform a biased random walk: if the bacterium senses it is moving up a concentration gradient towards an attractant it will rotate its flagella anti-clockwise more often than clockwise so as to continue moving in that direction; conversely, if it is moving away it is more likely to rotate its flagella clockwise so as to tumble with the aim of randomly changing its orientation to one which points it towards the attractant.

The chemotaxis of E. coli may be modelled (very) simplistically by considering a bacterium to move in a two-dimensional "world" populated by an attractant with a constant concentration gradient away from some location. At each of a series of time steps, a model bacterium detects whether it is moving up or down this gradient and either continues moving or tumbles according to some pair of probabilities.

Write a Python program to implement this simple model of chemotaxis for a world consisting of the unit square with an attractant at its centre. Plot the locations of 10 model bacteria which start off evenly spaced around the unit circle centred on the attractant location.

Solution P6.6.3

The code below tracks the bacteria positions over time on a two-dimensional plot (see Section 7.5.2 for a description of Matplotlib's imshow function.)

import numpy as np
import matplotlib.pyplot as plt

# The bacteria "world"
n = 100
world = np.zeros((n + 1, n + 1))

# Magnitude of step size
d = 0.02
# Number of bacteria, number of steps for each bacterium
nb, nsteps = 10, 400
# Position of the attractant
cx = cy = 0.5

# P0 is the probability of tumbling when heading toward the attractant; P1 is
# the probability of tumbling when heading away from it.
P0, P1 = 0.7, 1


def tumble():
    """Randomly pick a new direction, (dx, dy) for the bacterium."""

    theta = np.random.rand() * np.pi * 2
    return d * np.cos(theta), d * np.sin(theta)


def fpot(x, y):
    """Food potential, decreasing linearly as from (cx, cy)."""

    return -np.hypot((x - cx), (y - cy))


for b in range(nb):
    # Start the bacteria off evenly-spaced around the unit circle centred on
    # (cx, cy).
    phi = 2 * np.pi * b / nb
    x, y = cx + np.cos(phi) / 2, cy + np.sin(phi) / 2
    # Randomly pick an initial direction for the bacterium
    dx, dy = tumble()

    for i in range(nsteps):
        # We will estimate the gradient as proportional to the difference
        # between the food potential at the next position, (x',y') and that at
        # the current position: grad ∝ fpot(x',y') - fpot(x,y)
        grad = -fpot(x, y)

        # Update the bactrium's position (wrapping round at the world's edges
        x, y = max(min(x + dx, 1), 0), max(min(y + dy, 1), 0)
        ix, iy = int(x * n), int(y * n)
        world[ix, iy] = 1

        # So the attractant gradient at the new position is:
        grad += fpot(x, y)

        if grad < 0:
            # We're moving away from the attractant: definitely tumble
            prob = P1
        else:
            # We're not moving away, tumble only with this probability
            prob = P0
        # Decide whethe to tumble and update direction of movement if we do
        if np.random.rand() < prob:
            dx, dy = tumble()

plt.imshow(world, cmap="hot")
plt.show()

A typical plot is depicted below.