Team:XJTU-China/Expansion-Model
<!DOCTYPE html>
XJTU-China
Table of Contents
The Expansion Model
Background
Cellular Automaton (CA)
Microbial expansion in the environment has been an intriguing topic. However, it is too difficult to study accurate motion due to the complexity of a realistic system. Thus, theoretical models are urgently necessary.
The cellular automata model is proposed by Von Neumann in the 1950s to study the phenomenon of cell self-reproduction. It is a modeling method assuming that after setting the initial state, each cell is constantly updated according to certain evolution rules. CA is a system that is discrete in time, space, and state. Its basic elements include cellular, cellular space, cellular neighbor, and evolution rules.
Cellular
Cellular is the basic unit of cellular automata. A cell can be one-dimensional or two-dimensional, and its state is determined by a set of states at each discrete time point.
Cellular space
Cellular space is the space occupied by a distribution of cells. The appropriate shape can be selected according to the interaction relationship between cells in practical problems. The most common is the grid-type regular arrangement, including triangle, quadrilateral and hexagon.
Cellular neighbor
Cellular neighbors are a set of cells which have an effect on the change of cellular state. Generally, in the cellular automata model, the state of a cell is affected by the state of itself and surrounding cells, such as competitive nutrients or mutually beneficial symbiosis. In a numerical simulation, the range of neighbors can be determined by controlling the retrieval radius. Common types of neighbors are shown in Fig. 1.
Evolution rules
The evolution rules are set according to a practical problem. But generally, the state of the cell at time is a function of the state of the cell neighbor at time :
Where
- means time and is discrete;
- represents the state of the cell at time ;
- represents the state of neighbors of the cell ;
Our Model
Nishiyama et al.[1] applied the CA model into the study of the morphology of Bacillus Subtilis biofilm. Through simulation, they were able to observe the detailed process of B. subtilis expansion. To study our mutual system further, we want to propose a CA model based on their work, combining mutualism into the expansion process.
For the limit of time and current research, we could not gather all the effective parameters in the real world. Our goal is to:
implement a theoretical model that could be applied to a multi-species expansion system;
provide a method to predict the expansion of BSC;
observe the expansion intuitively and take a look at details.
Mathematical Model
Assumptions
To build up a CA model, the most important part is to set up evolution rules. Since we are simulating a growing, moving and interacting dynamic system, simplifications are necessary. Dividing the expansion into two phases each time step, we make the following principal assumptions:
Migration phase
The main forces driving microorganisms to move are[2]:
The osmotic pressure, which is formed by cell matter, EPS, other nutrient molecules, and ions. Areas with high osmotic pressure uptake water from the soil which causes swelling and expansion of the biofilm.
Cell “sliding” or “pushing”. When cells grow rapidly, they tend to occupy the whole space. As the population reaches $K_i$, they tend to push each other and the colony starts to expand outwards.
Factors we do not consider are:
- Cells don’t move themselves. Spores are inactive. Though B. Subtilis has flagella to move, they are downregulated when EPS is actively produced[2]. So the bacteria won’t remain motility that Nostoc sp. doesn’t have too.
We do not consider strong water flux. In the soil, water flux may take away EPS, cells, and nutrients [2]. Since there is little water in the sandy area, water flux is not strong enough to do this.
When updating states, we only consider the influence of a cell and its closest neighbors. In the later section, we calculate the “local gradient” by the difference of neighboring cells.
It’s natural because CPS is a component of the cell, while nutrients and RPS are not.
We assume that EPS (exopolysaccharide) is divided into RPS (released polysaccharide) and CPS (cellular polysaccharide)
Fig. 2 gives a clear illustration of the migration mechanisms.
Consolidation phase
Cell growth will be suppressed when they are dense. In our model, cells grow quickly and highly occupy the space (in order to “push”). It’s suitable to apply this basic assumption in Logistic growth. However, we assume that when cells grow actively (e.g when nutrients are enough), growth will not be totally suppressed.
About nutrients, we think
Nutrients refer to carbon sources in the soil. Nostoc sp. produces organic carbon while Bacillus Subtilis consume it.
Other necessary nutrients that Bacillus Subtilis provide for Nostoc sp. cannot be ignored. Without nitrogen and phosphorus, cyanobacteria could not survive. So the number of N & P solubilizing Bacillus Subtilis may decide whether Nostoc sp. can grow well.
Cells can be either vegetative cells or spores, and the two states can convert mutually. This is a property of both the microorganisms.
If nutrients are insufficient, both the microorganisms all change into spores instantly, while spores recover gradually only when there are abundant nutrients.
To specify terms, “cell” in later sections refers to a biological cell, while “cellular” means a parallelogram which contains multiple “cells”.
In our model, each cellular has 6 attributes:
| Attribute Number | Variable | Attribute Number | Variable |
|---|---|---|---|
| 0 | number of B. S | 1 | number of B. S spores |
| 2 | number of N. sp | 3 | number of N. sp spores |
| 4 | number of EPS | 5 | number of nutrition |
They are all nonnegative integers. Here we adopt a new system of units like atomic units:
unit nutrient: nutrient consumption of a B. Subtilis in a single timestep;
unit polysaccharide: the amount of CPS a B. Subtilis cell has.
unit length: the distance between the centers of a neighbor cell (in Fig. 1, unit length actually equals the side length of a cell in both Von Neumann and hexagonal grids)
We also have to specify several notations. Note: equals to "Attribute Number" listed in Fig. 1; means position index of a certain cellular.
| Symbols | Descriptions | remark |
|---|---|---|
| the number of the six attributes ( means time step) | ||
| the number of migration from idx to idx2 of the six attributes | ||
| the number of push-out from idx to idx2 of the two microorganisms | ||
| the number of the cells which turn into spores |
Here are the parameters, whose values are set by ourselves.
| Parameter | Descriptions | remark |
|---|---|---|
| birth rate | ||
| the probability to recover from spores | ||
| nutrition consumption per cell per time step | ||
| RPS production per cell per time step | ||
| carrying capacity | ||
| diffusion coefficient | ||
| amount of CPS a cell has | ||
| the number of B. S cells that can produce enough N and P for N. sp |
Evolution Process
Our evolution process is shown in this flow chart:
Initialize All Attributes
Put some bacteria, cyanobacteria, EPS and nutrition in some cells.
Calculate Gradient
We adopt Fick's first law, a classic equation to describe particles' diffusion in solution. Under the assumption that the diffusion flux () goes from regions of high solute concentration to regions of low concentration, with a magnitude proportional to concentration () gradient, we can write
in two dimensional space, it is
where is the diffusion coefficient, an intrinsic property of solute and the medium. A large value means higher diffusivity. is a 2D vector that can be decomposed into several directions. Therefore, we can calculate in any direction separately.
In this discrete model, we may use rather than the first derivative. is the difference of the number (of cells, EPS or nutrients) between neighbors. Under our unit system, . Thus, we can represent it as:
where is a direction vector of unit length.
Calculate Number of Migration
According to the quantitative relationship between osmotic pressure () and solute concentration () proposed by van't Hoff in 1886:
which means is proportional to . Since osmotic pressure has additivity, we can write:
According to assumption 1, cells, EPS and nutrition make up osmotic pressure with different proportional constants. We can compute "average osmotic pressure difference" by:
where is the weight. And we can formulate the number of migration into:
So for each attribution and each neighbor of a certain cell, we calculate the number of migration. Note that this value might also be negative.
Migrate and Update CPS
First, we update all attributions of all cells at the same time using . For each cellular,
According to Eq. 1~3, reverses if we exchange and . So the "income" term equals . We don't add the "outcome" term in order to avoid double counting.
Polysaccharide content is simultaneously updated according to the CPS assumption.
Update spores, nutrient and RPS
We first judge if the two microorganisms will generate spores. For B. Subtilis, if
which suggests carbon is insufficient, and for N. sp, if (recall that is a constant)
which suggests there is not enough B. Subtilis to produce N and P, then redundant cells turn into spores. For example,
In contrast, when nutrients are abundant, spores return to vegetive cells. For instance, when
spores of B. Subtilis recovers
We update the number of cells and spores: ( for all below)
where means the corresponding number of spores.
Then B. Subtilis consumes carbon nutrient while N. sp produces organic carbon:
Note that is negative. Finally, they produce RPS:
Cell Growth
We just use the Logistic growth model:
where is a space-limiting factor. When , the number of cells will decrease. However, we may multiply with a factor (slightly greater than 1) to release the growth limit.
Cell Pushing
After growth, the number of cells might exceed the carrying capacity. Redundant cells go out to its neighbor. They will certainly pick a place with fewer cells, in other words, more living space. So we assign weights for neighbors to receive redundant cells. The weight is also proportional to . So we get the number of cells to be pushed out:
where .
Then update cell states in random order.
Stop Criteria
When the simulation epoch reaches our preset maximum, or attribution values become negative, our program exits. Otherwise, the cycle goes on. We assign
In the whole program, we round non-integer attribution values after each step.
Result Analysis
We will show the results using a rectangular grid and "Moore" neighborhood. In a lack of realistic parameters, we set them by ourselves and try to give an intuitive explanation and fancy visualization of our model. The color bar gives us a relationship between color and
probs_migration = [0.02, 0.002, 0.002, 0.02] # diffusion coefficient
weight = [[1, 0.05, 1], [1, 0.05, 1]] # weight for osmotic pressure
params_BS = (0.40, 0.057, 0.1, 0.1, 1000) # growth rate, turn into spores, nutrient consumption, rps production, and carrying capacity
params_No = (0.01, 0.001, -0.1, 0.5, 1000) # N. sp grow much slower
init_values = (100, 200, 200, 1000) # B. S, N. sp, EPS, nutrient placed in the center
Fig. 4 shows the most ideal situation: enough nutrient (relatively low consumption) and a high growth rate. The colony forms a stable circle. At first, bacteria try to expand but fail for the nutrient limit. Then we can see a clear gradient (from yellow to dark green to purple) at an earlier time. Due to more and more rapid growth, grids in the center all reach the maximum carrying capacity, showing brilliant yellow. Only grids on the edge of the circle show a lower number of bacteria.
Fig. 5 shows how the distribution of EPS changes with time. Because EPS will not disappear, the highest value appears just in the center.
However, N. sp don’t grow due to a low growth rate. Thus, we set as 0.3 and run again. Fig. 6 shows that increasing the birth rate of Nostoc. sp results in a brilliant growth without negative effects on B. subtilis.
Under these parameters, we put the microorganisms (same amount) in four separate points as multiple colonies. As can be seen in Fig. 7, it doesn't affect every single colony's growth. With a tighter growth limitation (lower , the factor multiplied to ), the cells simply gather into a bigger circle, expanding outwards as a whole.
Fig. 8 shows the nutrient dynamics. At first, the nutrient is consumed quickly by the outer cells which grow actively. Everywhere they go suffer a nutrient decrease. But as more cells go through, nutrients can accumulate due to more production. After that, they are consumed by another group of expanding cells from the center.
The evolution of the number of spores is quite like the nutrient’s. The number of vegetative cells goes down on the edge first and is later made up by later microorganisms. We surprisedly find that spores might recover a lot with abundant nutrients. And there seems to be a cycle from more spores and less nutrient to fewer spores and abundant nutrient.
We also try to put more nutrients on the upper left of the figure to observe how the bacteria grow. We increase nutrient demand () to highlight the effect of nutrient limitation and run the program again. It seems that bacteria on the upper left corner has a little higher expansion rate, but there is no distinct difference. We guess it is because the two microorganisms provide enough nutrients for each other to grow and expand.
Conclusion and Prospects
We choose the CA model to describe the diffusion mechanism of a symbiotic system qualitatively. By introducing symbiotic mechanisms, we can further study the relationship between the two organisms.
The actual growth simulation of Bacillus subtilis is more suitable to use a hexagonal grid because in any direction the distance between two cells is the same (think of squares).
According to the literature reports, randomness might be important in morphology studies. And detailed environmental factors are considered. These are where we need to improve in the future.