"

From 2014.igem.org

Form a growing snowflake with Cellular Automaton

A growing snowflake can be programmed with two-dimensional cellular automaton[]. The cellular space is defined as a two - dimensional square lattice of cells with size of n*n. Every cell has two states, on and off, represented by “1” or “0” respectively, and has four orthogonally adjacent cells as neighborhood. Its evolution rule directs how cells change states from one discrete time point to the next and it is that: if a cell is off, it can change state to be on at the next time point only when one out of its four neighborhood is on; if a cell is on, it will remain on forever.

For a cell with coordination of , at time point , its state can be shown as. Then the evolution rule is:

Let initial state as all cells off except for the center four “on” cells. As time increases, a growing snowflake forms with less than 10 steps (one step equals to one time point). Figure 1 shows how a snowflake grows up at time point 6. Figure 2 shows a dynamic pattern of it.

Figure 1

Figure 2

In this experiment, the reaction process can be divided into three parts: signal receiving and signal feedback and signal transmission. The spreading of signaling molecules in culture satisfies certain regularity. The purpose of this model is to describe the molecular diffusion process, the process of signal transmission. Fick's law is one of the models described in physics molecular thermal motion, based on the fick's law to build the molecular diffusion model. This model is based on the following assumptions.

Hypothesis 1

No intermolecular interaction. Each molecule movement is independent of each other. The interactions between molecules is simplified in this model, so molecules move randomly in the medium, and the average moving rate is constant. So all the spread of the signal source is only related to the initial state and time, thus the problem is simplified to single source.

Hypothesis 2

There are two kinds of motions for molecular including vibration and spin at some point. The vibration frequency is f. The movement direction is random and it moves at an average rate and spins the original position unchanged.

Hypothesis 3

Molecular motion is always located on the same horizontal plane. It is convenient to make the movement rate in the horizontal direction to stand for the average speed. And it doesn't reflect molecules’ movement in the vertical direction.

This model uses one-dimensional cellular automata. The signal source is located in the first cell, and the cell number i stands for the distance between molecules and the signal source. Cellular signal source stands for circle medium which is in a i distance. For each of the molecules in a cell, there are three motion state at a same time, spin, i - 1 direction vibration and i + 1 direction, and the probability of vibration in different directions is equal. The signal molecules from signal source follow the secretion function of P, an average rate v and molecular vibration frequency f. This model is built depending on all above.

Optimization and improvement for the model is in the following.

There is interaction between molecules and the movement probability in the direction of I - 1 and I + 1 is influenced by the concentration. It is hard to measure the average rate of molecular motion and vibration frequency. And it may be better to use diffusion coefficient D instead.

Let’s view our project in a mathematical way.

To begin with, we need to do some simplification and give some assumptions for the system.
A) NO environmental factors will give any interference and stimuli to the system.
B) The whole culture plate for E.coli is perfectly in a shape of square and every E.coli colon occupies perfectly a square grid accordingly, accomplishing strictly square grids of lattice.
C) For every colon of E.coli, its depth is so equally minor that the diffusion of any molecular is restricted to two dimention.
D) Hence, the plate can be defined as a bi-dimentional space with orthogonal cartesian coordinate system. When we set the original point at the lower left corner of the vertices, and regulate the unit length as the side length of colon, we will use the coordinate of the higher right corner of the vertices as the colon’s coordinate.

1. Using system of ordinary differential equations to describe the genetic pathway

Among mathematical models for gene expression in the cell, ordinary differential equations played a dominant role in the history and are still widely used today despite the discovery of dramatic gene noise in the cell. Here, we use a system of ordinary differential equations to describe the interaction between the genetic parts in the designed genetic pathway of our project.

i) Assumptions:

According to the reference[] and fick’s law about diffusion process, we put forward the following assumptions for the model.
a) During the period of interest and under the specific environmental conditions, E.coli in the experiment is in good shape and remains steady expression of genes without proliferation;
b) No gene noise influences the function of the pathway;
c) The expression of genes in the pathway fits well in sigmoidal function ;
d) Both signaling molecules have no inter- and intra-molecular interaction and the diffusion of them obeys the Fick’s first law of diffusion, which postulates the diffusive flux is proportional to the concentration gradient;

ii) System of ordinary differential equations

Based on the assumptions above, we build two systems of ordinary differential equations, A and B, for the genetic pathway in the two species of E.coli respectively.
AHL：A;pC-HSL：B;CI:C;GFP：G;Complex of signaling molecule and LuxR or RpaR：R; concentration of LuxR：[LuxR]; Protein of Lac：L1,L2;environment :en;inside a cell:in;decay constant:γ;repressive constant:β.
To the species A E.coli, which can sense the signaling molecule AHL and produce the other kind of signaling molecule pc-HSL, the system A listed as below:

Similarly, to the species B E.coli, which can sense the signaling molecule pc-HSL and produce the other kind of signaling molecule AHL , the system B listed as below:

iii) Calculation and simulation in a discrete model

Considering we can not get the analytical solution for the system of ordinary differential equations above, we use a discrete model to calculate and simulate the whole system.
Our approach follows:
A) The discrete method functions as both spatially and temporally:
We firstly define the unit length of time as a relatively short time period , and it is understandable that the shorter the unit length is, the better the discrete model simulates the reality.
For a system with N*N colons, we divide each square colon into m*m minor grids of lattice, then there will be n*n (n=N*m) minor grids in the culture space. Then the unit length is defined as the side length of every minor square gird and equals to the diffusive distance during one unit length of time. (When we indicate that, we actually suggest that the diffusive velocity of one signaling molecule should be the integer times of that of the other.) An orthogonal coordination system sets up with the original point as the lower left corner of vertices as figure 1, and very minor grid gets a coordinate as . Every minor grid can be seen as homogenous.
B) To simulate that the two species of E.coli align alternatively in the culture plate, we define a discrete function system, determining which system of ordinary differential equations applying to the calculation for one minor grid .

C) For the concentration gradient at a grid , it is defined as:
D) Cellular automaton is used again here as a computation method for the system with iteration of multiple time units.

iv) Numerical simulation results

According to references, we set the value for parameters in the model; to make it even simpler, we assume the two quorum-sensing systems, AHL-LUXI/R in Vibrio fischeri pC-HSL-RPAI/R in Rhodopseudomonas palustris have identical functional and diffusive parameters. In this part, we use the production of GFP as reporter to demonstrate the performance of the designed genetic pathway and discuss the numerical simulation results.
A) A bell-shape curve for the GFP production to AHL concentration.

Figure 2

As figure 2, for any grid of E.coli colon, when put in different concentrations of targeting signaling molecule, they respond (means producing GFP) in a disparate efficiency and the responding curve is bell-shaped with a optimal concentration of signaling molecule.
B) The existence of a steady state for the production of both GFP and signaling molecules.

Figure 3-1

Figure 3-2

As the two figures shows, there exists a synchronized steady state for the production of both GFP and signaling molecule.
C)LacI decay constant, , positively correlating to GFP production.

Figure 4. The LacI decay constant, , under the 4 conditions is 0.2,0.4,0.6,0.8 respectively while other parameters are identical. The four curves come from the same steps of iteration.

As figure 4, the value of LacI decay constant, , positively relates to the production of GFP, which is rational because both of the two forms of LacI protein (with and without codon modified ) repress the Lac promoter, which control the production of GFP and signaling molecules. In addition, we can conclude that a bigger shortens the time length to the steady state.
D) The feasibility of low background expression of GFP and signaling molecules.

Figure 5. From down to up, steps of iteration of the system increase.

As figure 5 shows, after different periods of time, the background expression of GFP (corresponding to concentration of AHL equaling to 0) nearly remains the same. In addition, we can deduce a trend that whatever concentration of targeting signaling molecule is, a steady state will be achieved as long as enough steps of iteration performed.
E) The influence of concentration of LUXR to the optimal concentration of signaling molecules.

Figure 6. 3 peaks from left to right correspond to concentration of LUXR as 2,1,0.5 uM, respectively.

From figure 6, we can conclude that with the increase of the concentration of LUXR, slopes of the bell-shaped curves become steeper without the change of maximum expression of GFP during an identical period of time.
F) Pattern formation in multicellular system.
Use the MATLAB and C language, we have programmed and constructed this multicellular system. By modulating the parameters, initial state and iteration steps, we would accomplish multiple regular and dynamic patterns.
In figure 7-1, a complex dynamic pattern forms after 100 times of iteration with 1min as time unit. In figure 7-2, 7-3, another complex dynamic pattern forms with the increasing steps of iteration.
And in this pattern we modulate the parameters of the two species of E.coli to be slightly different and we get a pattern to be in symmetry but not central symmetry.

Figure 7-1. Planform (left) and side view (right) of a pattern. Two axis of [20,60] stand for the location coordination system of culture plate, and vertical axis of [0.2,0.5] stands for the production of GFP. The red color in the picture stands for relatively high concentration of GFP and the blue color in the picture stands for relatively low concentration.

Figure 7-2

Figure 7-3

This two figures show Planforms or side views of a pattern developing with the increasing steps of iteration of the system. The color bars beside demonstrate the relationship between concentration of GFP and color. In the development of this complex dynamic pattern, we can see a central square gradually develops through a complex rhombus stage , and towards two points on one diagonal line of the space.

Conclusion

In this model, we have constructed a multicellular system with complex dynamic patterns. By modulating its parameters, iteration steps and initial state , we can change its patterns. This model comes from a combination of systems of ordinary differential equations and cellular automaton and programmed with MATLAB and C language. Through the numerical simulation of program, we can use this model to predict and demonstrate complex multicellular system.