Modeling on Ecosystem Level

With the benefit of the DDEs model we built, we could acquire necessary data to move on to the Ecosystem level modeling based on Cellular Automaton. In this modeling, we intended to foresee the important points of the application. For example, we could determine how should we release Euk.Cement, how would Euk.Cement spread and how much Euk.Cement do we need to produce the best result. The ecosystem could derive the best strategy to make full use of Euk.Cement.

Part one: the “wake-up” problem:

Since it takes some time for the darkness induction system to shut down, we have to treat our strains in darkness to make sure they recover from the shut-down state. However, the problem is: How long should we put the strains into darkness?

As we can see from the DDEs model, our system has bi-stable state: the stable state in darkness and another stable state in light. The longer our strains are exposed to light, the longer it takes to convert them to the stable-state of darkness. Besides, our strains are more sensitive to light than darkness. In another words, if our strains are exposed to light for 10 minutes, it has to take 60 to 70 minutes for the strains to return to the original state.

For storage and transportation, we have to keep our strains exposed to light all the time to ensure they would not produce Si-tag or Mcfp-3. When we need them to work, we have to turn off the light and put them into darkness and wait for them to “wake up”. With the DDEs model, we can simulate the process of how do the strains “wake up”.

Results:

The figure shows that it takes about 3600 minutes for Euk.Cement to “wake up” from the stable-state of light. Therefore, before releasing our kit into the sands, we should turn off the light and treat them for 60 hours. If some Euk.Cement drift away from the sands, they will be exposed to sunlight and their downstream systems would be shut down by darkness induction system.

Part two: the permeation problem:

After releasing Euk.Cement, it would permeate the sands and we intended to find out what’s the final outcome of its permeation and what strategy should we take to attain the best result. To achieve this goal, we have built a permeation model with Lattice Method (Cellular Automaton)^[3] .

Hypothesis:

In order to make the model rigorous, we hypothesized the following prerequisites.
- There is no current in the sands.
- The total number of our yeast is constant. No yeast will be born or dead during the permeation process.
- The variability of temperature has no effect on our Euk.Cement or it could be ignored.
- The size of the sands and our Euk.Cement is constant.
- Every unit of sands has the same amount of silica.
The Cellular Automation simulated the process that our Euk.Cement expresses Si-tag and Mcfp-3 while it permeates in the sands. The PDEs (Partial Differential Equations) of this model are listed below.

Parameters:

The description of parameters, their values and the references involved in this model are listed in a table.

Parameter	Description
a(x,y,z,t)	the concentration of Mcfp-3
A(x,y,z,t)	the concentration of fixed bacterial
B	the vector of prerequisite coefficient
C(x,y,z,t)	the concentration of free bacterial
bi	the permeation coefficient of on the i direction
eb	the rate of Mcfp-3 expression
ed	the rate of Mcfp-3 degradation
l	the edge length of each lattice
T0	the time unit of prerequisite
α	the possibility of effective collision
λ	average length of steps of Brownian motion
τ	average free time of Brownian motion
NA	the Avogadro constant
T	the absolute temperature
R	the Rydberg constant
r & d	the radius and diameter of our yeast
ν	the speed of the prerequisite
ρ	the density of our yeast
η	the dynamic viscosity coefficient
Ci	the total number of bacterial of the 26 neighboring lattices of the lattice (x,y,z)

More detailed parameters table (2) can be downloaded here.

Formulary:

Monte Carlo Simulation:

To determine the values of some necessary parameters in the permeation model, we run the Monte Carlo simulation for a single lattice.
The Euk.Cement moves in the sands and has a certain possibility to bind the silica. Besides, the motion of Euk.Cement includes the Brownian motion, the sedimentation due to gravity and the motion due to the current. We hypothesized that there is no water flow in the sands, and therefore the motion of Euk.Cement is largely depended on the former two.

Brownian motion:

(2.6)

We know that the speed of Brownian motion: is relevant to the value of free time. The average length of steps of Brownian motion has the same scale with the cracks of the sands. Hence, we can set the unit time of Monte Carlo Simulation as 1s (It’s not the unit time of permeation). And with this setting, we could derive that every yeast has one collision with the sands.

Sedimentation due to Gravity:

(2.7)

and therefore

(2.8)

Comparing the effect of sedimentation due to gravity and the Brownian motion, we reach a conclusion that the latter one outweigh the former one (about 100 times). And therefore the most effective factor of the motion of Euk.Cement in the sands is Brownian motion.

The code of the Monte Carlo Simulation can be downloaded here. Its flow chart is

After running the Monte Carlo Simulation with 106 yeasts for multiple times, we finally determine the values of the following parameters

Lattices Method:

With sufficient parameters, we could finally simulate the permeation process. The code of our permeation model could be downloaded here.

We chose to use numerical solution to solve our Permeation model. We divided the space into the same lattices, and each lattice has only one phase. Hence, each variable in the same coordinate is considered the same. The timeline was divided into the minimum and discrete values. By calculating the variables in each lattice through the discretized timeline, we could attain the values of a(x,y,z,t),A(x,y,z,t)and C(x,y,z,t). If the lattices were small enough, the PDEs above could be replaced by Lattices Method. The equations of Lattices Method (Cellular Automaton) are listed below.

With the equations of Cellular Automaton as well as Boundary equations and Initial state equations, we could attain the values of the concentration of Mcfp-3, the number of fixed Euk.Cement and free Euk.Cement at any time and at any coordinate.

We ran the simulation with a server and printed a heat map

Results:

From the permeation process shown in the heat map, we can derive the following conclusion:
- The Euk.Cement permeates very efficiently over a large space. Therefore, perhaps we have to enclose it.
- A rock (1000 cm3) is added to the sands. However, it’s not even noticeable from the outcome. Hence, the rocks or stone inside the sands have little effect on the permeation.
- The final shape of the permeation is approximately global and not sensitive to the initial input.
- The Euk.Cement distributed quite evenly in the whole space they spread but not densely at the surface of the sands.
- The Euk.Cement expressed sufficient amount of Mcfp-3.

Robustness and Parameter Sensitivity Analysis:

The influential factors include noise of gene expression, variability of temperature, variability of the size of the sands and Euk.Cement, etc. However, our permeation model is most sensitive to the parameter α (the possibility of effective collision). Therefore, we focused on the analysis of the sensitivity of the parameter α.

When the α of MC = 0.001, the α of LM = 0.0282:
The total number of Euk.Cement bound with sands in 30s = 28224;
The total number of Euk.Cement drifted away in 30s = 74778;

val(;:;:1)=
	0	8	0
	4	12439	7
	0	4	0
val(;:;:2)=
	5	12474	4
	12456	1000000	12659
	2	12327	9
val(;:;:3)=
	0	7	0
	11	12423	4
	0	4	0

When the α of MC = 0.0013, the α of LM = 0.0365:
The total number of Euk.Cement bound with sands in 30s = 36464;
The total number of Euk.Cement drifted away in 30s = 74925;

val(;:;:1)=
	0	8	0
	3	12538	9
	0	7	0
val(;:;:2)=
	5	12553	3
	12549	1000000	12550
	2	12398	5
val(;:;:3)=
	0	5	0
	9	12337	4
	0	10	0

When the α of MC = 0.0015, the α of LM = 0.0421:
The total number of Euk.Cement bound with sands in 30s =42164;
The total number of Euk.Cement drifted away in 30s = 74661;

val(;:;:1)=
	0	2	0
	2	12619	2
	0	7	0
val(;:;:2)=
	7	12638	8
	12307	1000000	12394
	5	12351	6
val(;:;:3)=
	0	7	0
	8	12352	9
	0	4	1

Results:

From the result as shown above we could figure out that with 81.79% addition of the α, the bound Euk.Cement had an increase of 81.94%. This result indicates that the substance inside each lattice could be exchanged sufficiently to attain an even concentration of Euk.Cement as well as its secretions.
Therefore, we could change theα by selecting different Si-tag to achieve different targets, such as the intensity of the sands and the scope of the permeation. For example, we could change the scope Euk.Cement permeates and the intensity of sand by choosing different domain of Si-tag.

References:

[1] Chen, W.W.; Neipel, M.; Sorger, P.K. (2010). "Classic and contemporary approaches to modeling biochemical reactions". Genes Dev 24 (17): 1861–1875.doi:10.1101/gad.1945410. PMC 2932968. PMID 20810646.
[2] Kenneth Connors, Chemical Kinetics, 1990, VCH Publishers, pg. 14
[3] Wolf-Gladrow D A. Lattice-gas cellular automata and lattice Boltzmann models: An Introduction[M]. Springer Science & Business Media, 2000.
[4] Wolf-Gladrow D A. Lattice-gas cellular automata and lattice Boltzmann models: An Introduction[M]. Springer Science & Business Media, 2000.
[5] Sorokina O, Kapus A, Terecskei K, et al. A switchable light-input, light-output system modelled and constructed in yeast[J]. J Biol Eng, 2009, 3: 15.