I recently read an interesting review paper [Brown2009] that explores the possibility of fighting a bacterial infection with Darwinian medicine through the introduction of a “cheater” strain of bacteria into a wild population. I extended the analytic model in the paper to an agent-based model with explicit spatial structure.

This simulation incorporates three ideas from Section 2 of [Brown2009]:

• (a) A cheat that does not produce exoproducts (public goods),
• (b) Trojan horse cheats, and
• (d) Bacteriocinogen cheat invasion.

Instead of approximating spatial structure with the relatedness factor $r$ in 2(c) “Trojan horse cheat in a spatially structured host” I just explicitly put the bacteria in space – they diffuse in a two-dimensional continuum with periodic boundaries.

Combining (a), (b), and (d) gives the following dynamical equations for the rate of change of the wild ($W$) and cheat ($C$) strain densities: \begin{equation} \begin{array}{rcrccccl} \diff{W}{t} & = & W ( & \overbrace{1}^{\text{growth}} & \overbrace{- N}^{\text{competition}} & \overbrace{- x}^{\text{costs}} & \overbrace{+ b W / N}^{\text{public good}} & \overbrace{- a C / N}^{\text{toxication}} & \overbrace{- e C / N}^{\text{bacteriocide}} )\\ \diff{C}{t} & = & C ( & 1 & - N & - q & + b W / N & - a C / N ) \end{array} \label{eq:qssa} \end{equation} where

• $N = W + C$ is the total bacterial density,
• $x$ is the cost of producing the public good,
• $q$ is the direct growth cost of the cheat (if any),
• $b$ is the benefit of the public good,
• $a$ is the strength of the antibacterial toxin, and
• $e$ is the efficacy of the bacteriocinogen.

These two equations combine Eqs. (2.1)-(2.5) from [Brown2009] – the effects can be explored separately in the simulation by moving sliders of other parameters to zero. The numeric dynamics are displayed as faint curves in the “dynamics” and “phase” graphs. They are included to validate and compare with the agent-based model.

The numeric model can be derived as a non-spatial limiting case of an agent-based model. In this model we explicitly represent each individual “agent” in the system instead of treating them as identical and just tracking aggregate population densities. A convenient approach is to borrow reaction kinetics from physical chemistry to describe allowed interactions in our model.

To capture the numeric model we require five types/species¹⁰⁾:

• $W$ is a wild-type bacterium,
• $C$ is a cheat bacterium,
• $F$ is a parcel of food,
• $T$ is a parcel of toxin, and
• $B$ is a parcel of bacteriocinogen.

We simply allow these agents to interact stochastically at given rates according to the following reactions: \begin{equation} \begin{array}{rclrcll} N & \xrightarrow{1} & 2 N & & & & \text{(growth)} \\ W + N & \xrightarrow{1} & N, & C + N & \xrightarrow{1} & N & \text{(competition)} \\ W & \xrightarrow{x} & \emptyset, & C & \xrightarrow{q} & \emptyset & \text{(costs)} \\ W & \xrightarrow{b} & W + F, & N + F & \xrightarrow{\beta} & 2 N & \text{(public good)} \\ C & \xrightarrow{a} & C + T, & N + T & \xrightarrow{\delta} & \emptyset & \text{(toxication)} \\ C & \xrightarrow{e} & C + B, & C + B & \xrightarrow{\gamma} & C & \text{(bacteriocide)} \\ & & & W + B & \xrightarrow{\gamma} & \emptyset \end{array} \label{eq:agent} \end{equation} where $\emptyset$ indicates an absence of products and $N$ is shorthand for either $W$ or $C$.

A moment's review should satisfy the reader that these reactions are reasonable prescriptions. For example, $W \xrightarrow{b} W + F$ indicates that each wild-type bacterium produces a food parcel at an average rate of $b$ per unit time. A bacterium that encounters a food parcel consumes it at rate $\beta$ and benefits by reproducing, $N + F \xrightarrow{\beta} 2 N$.

In the agent-based approach we model these processes explicitly and track all of the agents as they are produced, changed, or removed from the system. To connect with the numeric model we apply the law of mass action to compute the dynamics of a large, well-mixed population following the above reactions: \begin{equation} \begin{array}{rcrccccl} \diff{W}{t} & = & W ( & \overbrace{1}^{\text{growth}} & \overbrace{- N}^{\text{competition}} & \overbrace{- x}^{\text{costs}} & \overbrace{+ \beta F}^{\text{public good}} & \overbrace{- \delta T}^{\text{toxication}} & \overbrace{- \gamma B}^{\text{bacteriocide}} )\\ \diff{C}{t} & = & C ( & 1 & - N & - q & + \beta F & - \delta T )\\ \diff{F}{t} & = & & & & & b W - \beta F N\\ \diff{T}{t} & = & & & & & & a C - \delta T N\\ \diff{B}{t} & = & & & & & & & e C - \gamma B N. \end{array} \label{eq:full} \end{equation}

Eq. ($\ref{eq:full}$) looks much more complicated than ($\ref{eq:qssa}$) we're trying to relate it to. But we can reduce the system if we assume that the densities $F$, $T$, and $B$ are fast variables – that they respond quickly to perturbations and quickly converge to equilibrium. If we assume that they converge so quickly that the slow variables ($W$ and $C$) don't change appreciably while the fast ones equilibrate then we can approximate the fast variables as always being in equilibrium – the quasi-steady-state approximation (QSSA): \[ \begin{array}{rcl} F & = & b W / \beta N\\ C & = & a C / \delta N\\ B & = & e C / \gamma N. \end{array} \] Then Eq. ($\ref{eq:full}$) reduces to exactly ($\ref{eq:qssa}$) and the agent-based model in Eq. ($\ref{eq:agent}$) is an extension of [Brown2009] to finite populations. Unfortunately, finding the conditions to satisfy QSSA is not trivial [Segel1989]. Nevertheless, it should be satisfied if all of the rate constants $\beta$, $\delta$, and $\gamma$ are large. In the simulation I set them all to the same large constant so we may expect differences between the agent-based and numeric models result from other causes, such as spatial structure.

Brown et al. [Brown2009] create a structured population in Section 2(c) through preferential interactions of the bacteria with kin, mimicking spatial segregation. Implementing it as an agent-based model in [NetLogo] allows me to explicitly include space: agents move continuously in a two-dimensional plane with periodic boundary conditions. Reactions are localized to patches on a square grid. In the limit of large population size and infinitesimal patch size mass action swamps stochastic effects and the model should become equivalent to a reaction-diffusion system.

I arbitrarily chose to allow only the bacteria (wild and cheater types) to move. Byproducts (food, toxin, & bacteriocinogen) are stationary after being created. Some other options would be to allow the byproducts to move (for example, as if they were diffusing via Brownian motion) or for both bacteria and byproducts to move (possibly at different rates).

Unlike other events, which are treated as probabilistic stochastic processes, every moving agent moves in every time increment according to a random walk, with a jump size is governed by the diffusion constant. The effects of spatial structure can be minimized by mixing¹¹⁾. Only bacteria move so there may still be some spatial effects from the stationary byproducts but the dynamics should be more similar to the numeric (mean field) approximation.

A Trojan horse approach to medical intervention strategies

A NetLogo model by Rik Blok.

http://www.zoology.ubc.ca/~rikblok/wiki/doku.php?id=science:popmod:brown2009:start

This agent-based model represents a bacterial infection as described in [Brown2009]. It is assumed the wild type (W) produces a public good (F) which benefits all. To combat the infection the population is innoculated with a "cheater" strain (C). Changes are simply represented as elementary reactions between (local) agents:

• N → 2 N @ rate 1 (growth)
• W + N → N @ rate 1 (competition) C + N → N @ rate 1
• W → W + F @ rate b (public good) W + F → 2 W @ rate β
• W → ∅ @ rate x (wild cost)
• C → ∅ @ rate q (cheater cost)

where N indicates an individual of either type and ∅ indicates an absence of products.

Additionally, the cheater strain has one or more of the following traits:

• It consumes but does not produce the public good (F), C + F → 2 C @ rate β
• It produces a toxin (T) that harms all, C → C + T @ rate a N + T → ∅ @ rate δ
• It produces a bacteriocinogen (B) it is immune to but harms the wild type, C → C + B @ rate e C + B → C @ rate γ W + B → ∅ @ rate γ.

The target and infected cells are fixed whereas the virions move randomly with diffusion constant diff-const. If 'well-mixed' is on then they are shuffled randomly in space with each timestep.

Finally, it is assumed that the wild type is more "virulent" (harmful to the host) than the cheater strain.

How it works

The simulation approximates a Poisson process for each of the above events. The best known technique would be the Gillespie algorithm [Gibson2000] but it isn't well suited to NetLogo's strengths. Instead, time proceeds in steps with multiple events occurring in each timestep.

The step size is adaptive, chosen to achieve a desired error tolerance, compared with the Gillespie algorithm. When the error tolerance is near zero the likelihood of each event is small and we may expect just a few events to occur per timestep. Then we're accurately -- but inefficiently -- mimicking the Gillespie algorithm. As the tolerance increases we have more simultaneous events, lowering accuracy but increasing performance.

References

[Brown2009] Brown, Sam P., West, Stuart A., Diggle, Stephen P., and Griffin, Ashleigh S. Social evolution in micro-organisms and a Trojan horse approach to medical intervention strategies. Philosophical Transactions of the Royal Society B: Biological Sciences, 364: 3157-3168. doi:10.1098/rstb.2009.0055. 2009.

[Gibson2000] Gibson, Michael A. and Bruck, Jehoshua. Efficient Exact Stochastic Simulation of Chemical Systems with Many Species and Many Channels. J. Phys. Chem. A, 104(9): 1876-1889. doi:10.1021/jp993732q. 2000.

[NetLogo] Wilensky, U. NetLogo. http://ccl.northwestern.edu/netlogo/. Center for Connected Learning and Computer-Based Modeling, Northwestern University. Evanston, IL. 1999.

Like this page? Share it!

• Twitter
• Facebook
• Google+
• LinkedIn
• Pinterest
• Tumblr
• Reddit
• Email