A coarse-grained model of single particle diffusion
Part of a modeler’s job is to find simple models that capture the essence of a system while running quickly and scaling well to larger inputs.
In our case, we have a very “fine-grained” reaction-diffusion model illustrating Turing patterns, but this model consumes a huge amount of computational resources because it requires tracking the movements of hundreds of thousands of individual particles. Our goal is to build a “coarse-grained” model that will allow us to appreciate Turing patterns without the computational overhead required to track particles.
Our idea is to grid off two-dimensional space into blocks and store only the concentration of each type of particle found inside the block. To simplify things further, we assume that there is some maximum concentration of particles possible, so that we can divide the number of particles by this maximum concentration. As a result, the concentration of a particle in each block will be represented by a decimal number between 0 and 1.
Let us begin with an example of the diffusion of only A particles; we will later add B particles as well as reactions to our model. Say that the particles are at maximum concentration in the central cell of our grid and are present nowhere else, as shown below.
A 5 x 5 grid showing hypothetical initial concentrations of A particles. Cells are labeled by numbers between 0 and 1 representing their concentration of A particles. The central cell has maximum concentration, and no particles are contained in any other cell.
We will now update the grid of cells after one time step to mimic particle diffusion. To do so, we will spread out the concentration of particles in each square to its eight neighbors. For example, we could assume that 20% of the current cell’s concentration diffuses to each of its four adjacent neighbors, and that 5% of the cell’s concentration diffuses to its four diagonal neighbors. Because the central square in our ongoing example is the only cell with any particles, the updated concentrations after a single time step are shown in the following figure.
Note: The sum of the values in the grid in the figure above is 1, which is the same as the sum of the values in our original figure. Regardless of how many times we update the grid, the sum of the values should be 1 to ensure the conservation of total mass in the system.
After an additional time step, the particles continue to diffuse outward. For example, each diagonal neighbor of the central cell in the above figure, which has a concentration of 0.05 after one time step, will lose all of its particles in the following step. Each of these cells will also gain 20% of the particles from two of its adjacent neighbors, along with 5% of the particles from the central square (which doesn’t have any particles). This makes the updated concentration of this cell after two time steps equal to 2(0.2)(0.2) + 0.05(0) = 0.08.
The four cells directly adjacent to the central square, which have a concentration of 0.2 after one time step, will also gain particles from their neighbors. Each such cell will receive 20% of the particles from two of its adjacent neighbors and 5% of the particles from two of its diagonal neighbors, which have a concentration of 0.2. Therefore, the updated concentration of each of these cells after two time steps is 2(0.2)(0.05) + 2(0.05)(0.2) = 0.02 + 0.02 = 0.04.
Finally, the central square has no particles after one step, but it will receive 20% of the particles from each of its four adjacent neighbors, as well as 5% of the particles from each of its four diagonal neighbors. As a result, the central square’s concentration after two time steps is 4(0.2)(0.2) + 4(0.05)(0.05) = 0.16 + 0.01 = 0.17.
In summary, the central nine squares after two time steps are as shown in the following figure.
STOP: What should the values of the “?” cells be in the above figure? Note that these cells are neighbors of cells with positive concentrations after one time step, and so their concentrations should be positive after two time steps.
The coarse-grained model of particle diffusion that we have built is a variant of a cellular automaton, or a grid of cells in which we use fixed rules to update the status of a cell based on its current status and those of its neighbors. Cellular automata form a rich area of research applied to a wide variety of fields dating back to the middle of the 20th Century; if you are interested in learning more about them from the perspective of programming, then you might like to check out the Programming for Lovers project.
Slowing down the diffusion rate
There is just one problem: our diffusion model is too volatile! The figure below shows the initial automaton as well as its updates after each of two time steps. In a true diffusion process, all of the particles would not rush out of the central square in a single time step, only for some of them to return in the next step.
Our solution is to slow down the diffusion process by adding a parameter dA between 0 and 1 that represents the rate of diffusion of A particles. Instead of moving a cell’s entire concentration of particles to its neighbors in a single time step, we move only a fraction dA of them.
Revisiting our original example, say that dA is equal to 0.2. After the first time step, only 20% of the central cell’s particles will be spread to its neighbors. Of these particles, 5% will be spread to diagonal neighbors, and 20% will be spread to adjacent neighbors. The figure below illustrates that after one time step, the central square has concentration 0.8, its adjacent neighbors have concentration 0.2dA = 0.04, and its diagonal neighbors have concentration 0.05dA = 0.01.
Adding a second particle to our diffusion simulation
We now will add B particles to the simulation, which we assume also start with 100% concentration in the central square. Recall that B, our “predator” molecule, diffuses half as fast as A, the “prey” molecule. If we set the diffusion rate dB equal to 0.1, then our cells after a time step will be updated as shown in the figure below. This figure represents the concentration of the two particles in each cell as an ordered pair ([A], [B]).
A figure showing cellular concentrations after one time step for two particles A and B that start at maximum concentration in the central square and diffuse at rates dA = 0.2 and dB = 0.1. Each cell is labeled by the ordered pair ([A], [B]).
STOP: Update the cells in the above figure after another generation of diffusion. Use the diffusion rates dA = 0.2 and dB = 0.1.
Visualizing particle concentrations in an automaton
As we move toward diffusing a large board that is hundreds of cells wide, listing the concentrations of our two particles in each cell will be difficult to analyze. Instead, we need some way to visualize the results of our diffusion simulation.
First, we will consolidate the information stored in a cell about the concentrations of two particles into a single value. In particular, let a cell’s particle concentrations be denoted [A] and [B]. Then the single value [B]/([A] + [B]) is the ratio of the concentration of B particles to the total number of particles in the cell. This value ranges between 0 (no B particles present) and 1 (only B particles present).
STOP: What should be the value of [B]/([A] + [B]) if both [A] and [B] are equal to zero?
Next, we color each cell according to its value of [B]/([A] + [B]) using a color spectrum like those shown in the figure below. We will use the
Spectral color map, meaning that if a cell has a value close to 0 (relatively few predators), then it will be colored red, while if it has a value close to 1 (relatively many predators), then it will be colored dark blue.
When we color each cell over many time steps, we can animate the automaton to see how it changes over time. We are now ready for the following tutorial, in which we will implement and visualize our diffusion automaton using a Jupyter notebook.
The video below shows an animation of a 101 x 101 board with dA = 0.5 and dB = 0.25 that begins with [A] = 1 for all cells. All cells have [B] = 0 except for an 11 x 11 square in the middle of the grid, where [B] = 1. Without looking at individual concentration values, this animation allows us to see immediately that the A particles are remaining in the corners, while a band of B particles expands outward from the center.