## Overview

We cannot directly measure the number of workers, cheaters and guards in the mixed condition in the laboratory，because it is difficult to distinguish them. To verify the effectiveness and stability of cheat me if you can in industrial or other application scenarios, we choose to use mathematical tools to model the growth of the three bacteria.

Our goals:

1. Use mathematical models to simulate the dynamic behavior of cheaters, guards and workers competing with each other under industrial conditions, verify the effectiveness and stability of the platform.

2. Calculate the optimal value of the guard input under industrial conditions to maximize the output of the target substance.

3. Deepen our understanding of the suppression mechanism of guards on cheaters.

## Modeling Synthesis

To model the output of target substances under the condition of using our platform, we simulated the competition of cheaters, guards and workers. Our main model is the logistic competition equation, and in which we designed three models to simulate the killing of the cheater by the guard: the cheaters produce AIP, the AIP stimulates the guard to produce mCherry, and mCherry induces the cheater to commit suicide. We use the above process to model the killing effect of the input of the guard on the cheater, and obtain the killing rate of the cheater in the competition equation.

## Basic Model Assumptions

###### 1. The total volume of the system remains unchanged.

Our experimental environment is solution. If the total volume of the system is assumed to be constant, the concentration can be directly used to express the amount of substances and bacteria, ignoring a constant multiple.

###### 2. The concentration of public substances, toxins, and densities of different bacteria in different positions in the solution are even.

Due to constant stirring in practice, we made this assumption. This assumption eliminates the space factor and simplifies the model.

###### 3. Nutrients are single in our system.

Under this hypothetical condition, the competition among worker, cheater and guard is due to their consumption of common substances. Therefore, we can evaluate the competitiveness through the consumption of sugar in the experiment.

###### 4. The cheater mutates from the worker, and the mutation probability is constant.(the value is between 10-8 and 10-3)

###### 5. The concentration of the target substance is proportional to the concentration of the worker.

The quantity of target substances produced by the guard in the system is so small that it can be ignored, only the workers are considered.

###### 6. The guard’s ability to secrete target substances and mCherry does not affect its growth rate and environmental capacity.

When the guard produces mCherry, its pathway to synthesize the target substance is inhibited. Therefore, we approximately ignore the changes in the metabolic pressure of the guard.

###### 7. Ignore the process of substances entering and leaving the cell.

It is difficult to distinguish the fluorescence intensity inside and outside the cell in real time during the experiment, so this assumption is to match the experimental data.

###### 8. The production rate of AIP is proportional to the concentration of the cheater.

AIP is only 7 amino acids long, and we think its synthesis speed is very fast.

###### 9. Both AIP and mCherry will degrade naturally.

We found that if the initial concentration of AIP is given, the fluorescence intensity of mCherry rises first and then falls within one hour through experiments, proving that mCherry will degrade. When the initial value of AIP increases in a very large range (5nM-25000nM), the time when the mCherry fluorescence intensity reaches the peak is not significantly delayed, so we can expect that AIP may naturally degrade. The fitting effect after considering the degradation of AIP is better than not considering it.

###### 10. The concentration of AgrA produced by the guard is adequate.

The process of phosphorylation of AgrA by AgrC* into AgrAP can be described by the Michaelis-Menten equation, but it is difficult for us to obtain experimental data on concentration changes in this process. Based on this assumption, we get that the rate of AgrAP production is proportional to AgrC*, which simplifies the equation and reduces the parameters to be fitted.

###### 11. The kill rate of mCherry to the cheater is equal to the binding rate of mCherry to the surface receptor of the cheater.

The cheater will definitely be killed as soon as the suicide protein is activated. It is just a matter of time. The surface receptors of the killed cheater will be inactivated and cannot bind to mCherry. Therefore, we assume that the kill rate of mCherry to the cheater is equal to the binding rate of mCherry to the surface receptors of the cheater.

## Without platform

### The competition model

In the absence of a platform, the worker(W) and the cheater(C) conform to the logistic equation when they live alone in a natural environment. Considering that the worker will mutate into the cheater with a certain probability, we can get the following equation:^{[1]}

Here, the parameter p is given by us, the growth rate r and environmental capacity K is derived from fitting experimental data with logistic equation. And glucose assumption ratio is directly from wet lab.

Fig.1. The fitted growth curve

We can see in Fig2 that even with a very small mutation probability p, Cheater will totally replace Worker in about 70h.

Fig. 2. The dynamic competition behavior between workers and cheaters without platform (using OD600 to measure the concentration of bacteria)

## With platform

### 1.The competition model

Although the r and K change in the presence of the platform, the Cheater will still achieve an overwhelming victory over the competition with the Worker, as you can see in Fig.3

Fig. 3. The simulation of Cheaters and Workers in the presence of the platform.

Under the condition of introducing the platform, we added guards in the system. When the worker(W), cheater(C), and guard(G) are living alone in a natural environment, we assume that the growth of their quantity follows the logistic equation. Guards can not only compete with Workers and Cheaters, but also kill Cheaters, which is denoted by variable θ.

The equation we used to model the relationship between the killing rate and the concentration of mcherry is Hill equation in pharmacodynamics.^{[2]}

### 2.The Kill Model (Estimate of θ)

When there is the cheater in the solution, it will produce AIP, and AIP will stimulate the guard to produce mCherry. Then mCherry will bind to the surface receptor of the cheater to induce the cheater to commit suicide. The rate at which the cheater is killed is related to the concentration of mCherry in the solution.

To describe the above process, we will model the killed rate of the cheater in three steps.

###### Step 1: The cheater product AIP

In order to estimate AIP production rate, we communicated with the experimental team and did the following two experiments.

Firstly, we put different concentrations of AIP into the bacterial solution of the guard, and then measure the amount of mCherry produced by the guard. The above data models the relationship between AIP and the peak value of mcherry produced by the guard.

Fig.4. The relationship between the Fluorescence intensity of GFP and AIP(mmol/L)

We found that the two have a very good linear relationship, the expression is y=0.0011x+13.667, and R^{2 is 0.8055. Therefore, the initial AIP concentration that stimulates the guard can be approximated by the amount of mCherry produced by the guard.
}

Second, we let the cheater generate 10min and 15min of AIP, and then put the AIP into the bacteria solution of the guard to measure the amount of mcherry produced. Substituting the amount of mCherry into the above expression, we can deduce the amount of AIP corresponding to the two time nodes. By subtracting the two, we can deduce the amount of AIP produced by the spoofer within 5 minutes, and then get the rate of AIP produced by the cheater. The cheater per OD produces 1,369.86 nM of AIP per minute.

###### Step 2: The AIP stimulates the guard to produce mcherry

Fig.5. The path diagram of kill model

The free AIP in the solution will bind to the receptor AgrC on the surface of the guard cell, and activate AgrC into AgrC*. When Agrc is activated to AgrC*, AgrC* acts like an activating enzyme to phosphorylate AgrA inside the guard into AgrAP. When AgrAP binds to the promoter, mRNA will be transcribed, mRNA will be translated into mCherry, and mCherry will be secreted to the cell. Outside.The changes of substances in this process are as follows:

Ignoring the time it takes for AIP to bind to Agrc receptors on the surface of guard cells, assuming instantaneous equilibrium, the Hill equation can be used to model this process. We also need to consider the natural degradation of AIP (this is found to be necessary during fitting). The specific expression is shown below:

According to the basic assumptions mentioned above, the production rate of AgrAP is proportional to AgrC*:

In order to achieve a better fitting effect, the transcription and translation equations of mCherry are expressed separately. Considering the hydrolysis of mRNA and mCherry, the equations are as follows:

We estimate the parameters in Equation from experimental data: Given the amount of guard, it is believed that it will not change over time. Put different initial concentration gradients of AIP into the solution, and measure the change of mCherry's fluorescence intensity over time. In this case, since there is no cheater to produce AIP, the AIP in the system will only decrease as it degrades after being combined with Agrc, so the formula (1) is changed as follows:

Other formulas remain unchanged.

With known initial values of [AIP],for any given set of parameters the variation of [mCherry] with time can be solved by equations (2)-(6) (for corresponding to the experimental time conveniencely, we have used here an explicit Eulerian format with a time step of 0.01 min),and the error s is given by

We use a genetic algorithm to estimate the parameters in equations. Each individual in the population represents a set of parameter values, the corresponding s is the fitness of the individual, and finally the parameter value that minimizes s is selected.

The results of parameter estimation are as follows:

The fitting result of the theoretical curve and the experimental data is shown in the figure 6.

Fig.6. Fitted curves(with genetic algorithm)

Note that n in the fitting result is greater than 1, which indicates that the combination of AIP and guard has a positive synergistic effect.

###### Step 3: mcherry induces the cheater to commit suicide

In order to measure the killing ability of mCherry generated by the guard against the cheater, we designed a new pathway. Originally, mCherry would bind to the receptor on the surface of the cheater and activate the suicide gene ccdb of the cheater, thereby inducing the cheater to produce suicide protein and commit suicide. To facilitate detection, we replaced the ccdb gene with the green fluorescent protein(GFP) gene.

Fig .7. a new pathway

The cheater will definitely be killed as soon as the suicide protein is activated. It is just a matter of time. The surface receptors of the killed cheater will be inactivated and cannot bind to mcherry. Therefore, we assume that the kill rate of mCherry to the cheater is equal to the binding rate of mCherry to the surface receptors of the cheater.

In addition, the number of GFP produced is approximately proportional to the number of receptors bound by mCherry. Therefore, we put different concentrations of mCherry into the cheater with a fixed concentration, and then measured the GFP fluorescence intensity produced by the cheater, and used the fluorescence intensity to express the killing intensity of the cheater.

Since the value of GFP is basically stable after 10 minutes, we take the GFP fluorescence intensity corresponding to the experiment of different mCherry concentrations at 10 minutes to express the amount of GFP.

The equation we used to model the relationship between the killing rate and the concentration of mcherry is Hill equation in pharmacodynamics to model the relationship between mCherry and kill rate.

Among them θ is the fraction of receptor protein that has been occupied by the ligand, that is, the receptor binding rate. L is the concentration of free (unbound) ligand, that is, the concentration of mCherry. k_d is the equilibrium constant for dissociation. And n is the Hill coefficient. In order to facilitate the fitting, we simplify the equation, take the logarithm, and integrate, and finally get the following equation:

First, we use GFP to express the binding rate of the receptor bound by mCherry.First, standardize GFP to make it between [0,1], and the changes are as follows:

Among them, c and k are parameters, c represents the kill rate corresponding to the minimum value, and (c+1/k) represents the maximum kill rate.Then substitute the integrated Hill equation for fitting.When c and k change, the residual error corresponding to the experimental data is shown in the figure:

Fig.8. Residuals when taking different parameter values

We find that when the interval is [0.0011,0.3333], the fitting effect is the best, that is, when c=0.0011, k=3, the fitting effect is the best. At this time, the Hill equation parameters are: n=2.2347, kd=exp (10.2764).

The corresponding Hill equation image is as follows:

Fig.9. Hill equation fitting of GFP with mcherry

It is noted that the Hill coefficient n is greater than 1, indicating that the binding of mCherry and cheater has a positive synergistic effect. And in the communication with the experimental group, we got that mCherry can be observed to have clumping phenomenon by fluorescence microscopy during the experiment, which may be an explanation for the fitting effect.

### 3.Results

By numerically solving these equations, we found that because of the strong competitiveness, the Guard will not only eliminate Cheater, but also inhibit the growth of Worker, see in Fig10.

So how can we make this platform work?

Fig. 10. The simulation of platform’s effect when Guards are mixed with Cheaters and Workers in the system. The guard will not only eliminate cheaters, but also inhibit the growth of workers.

## Consider the condition of Alginate beads

By communicating with the wet lab, we decided to use Alginate beads to artificially limiting Guards, so that both the elimination of Cheaters and the domination of Workers can be guaranteed.

### 1. Modified model

It is assumed that the function of the Alginate beads is to reduce the guard's environmental capacity in proportion to the volume through space constraints, thereby intensifying the internal competition of the guards. However, the existence of Alginate beads does not affect the diffusion of substances, that is, the materials inside and outside the Alginate beads are evenly distributed, and the guards still compete with the cheaters and workers by consuming glucose, so it does not affect the competition between the guards, workers and cheaters. In a word, the guard's ability to compete with workers and cheaters remains unchanged, but its own environmental capacity will decrease.

Under the assumption that the ratio of the volume of the Alginate beads to the total volume of the system is α, the equations of the worker and the cheater remain unchanged, and the guard competition equation is changed as follows:

### 2. Results

You can see the simulation here in Fig11.

Fig. 11. The simulation of platform’s effect with Sepharose beads, where p=0.001, a=0.054

After introducing the platform, we found that the amount of worker would be maintained at a higher level( nearly environmental capacity), while the amount of cheater would be maintained at a lower level, which is lower than 2% of worker. indicating that our platform could indeed suppress cheater and increase the production of target substances. The simulation verified the effectiveness of our platform.

### 3. The optimal guard ratio

It’s obvious that for a given ‘p’, the final station of this system will change with varying ‘a’, see in Fig 12, where the best ‘a’ which guarantees the maximum of Workers is about 5%,that’s the value we need.

Fig.12. Final stations of the system with varying values of the input guards , where p=0.001

Notice that when ‘a’ is larger than ‘best a’, the amount of Worker will lightly decline. However, when ‘a’ is smaller than ‘best a’, the amount of Worker will decline sharply.

Therefore, the best value of parameter ‘a’ with varying mutation probability ‘p’ is shown in table3 and Fig13.

Fig.13. The optimal Guard ratio‘a’ with varying mutation probability ‘p’.

We found that when the mutation probability is less than or equal to 1 * 10 -3 , the amount of guard input needs to be increased from about 1% to 5% as the mutation probability increases. And the amount of worker slightly declines by 15%, but is relatively stable. When the mutation rate is 1 * 10 -2 , the guard will not be able to suppress the cheater, the platform will no longer work, and the amount of worker will eventually tend to 0. See in Fig.14.

Fig.14. The platform will no longer work

This gives the range of applicability of our platform and the optimal amount of guard input needs corresponding to different mutation probabilities.

## Sensitivity analysis

We analyze the sensitivity of our results by changing the value of parameters.

### The target gene

We found that the result is robust to the growth rate “r” and the environmental capacity “K”，which also influence the competition item. It means that when we change the target gene, our platform is still effective.

Fig.15. The result is robust to the growth rate “r” and the environmental capacity “K”.

When the growth rates was changed, the change of the results is not of significance.f or a disturbance with 10 percent, the maximum fluctuation of the Worker end point is only 4.7 percent, while that of the optimal guard input is only 1 percent.

Now let’s look at a disturbance in the environmental capacity. if the environmental capacity of the Guard is varying 10 percent, there is no significant change in optimal guard input. However, there will be a relatively huge fluctuation for the workers’ end point. That is probably because of the growth rate of The guard is big that a small changed in the environmental capacity will cause a strong competition between workers and guards, causing change in workers’ end point. That is the same reason for the big fluctuation in the optimal Guards input while a disturbance in workers’ environmental capacity. The model is robust for Cheater’s environmental capacity.

### The affinity of receptors to mCherry

However, when the hill constant of the combination of mCherry and cheater changes, the result of the model changes greatly. Specifically, when n becomes larger, the change is small. But when n decreases, as shown in the figure, there is a great change from the previous results. It’s importance to notice that hill constant represents the binding ability of mCherry and guard according to documents. When n > 1, it means that the binding is positive synergy, while when n < 1, the binding is negative synergy. Our fitting n is 2.247. When n decreases, that is, when n approaches 1, it means that the binding mode may change, which is naturally quite different from our original results. However, if only for the case of increasing n, the model is still robust.

Furthermore，we noticed that if we enhance the affinity of receptors to mCherry ( by decreasing Kd) ， the best amount of guard input will be significantly reduced，which means our platform will be more efficient.

Fig.16. if we enhance the affinity of receptors to mCherry， our platform will be more efficient.

The relationship here is roughly exponential, as you can see in the figure.

## References

[1]Zhao Chunxi. Three Group Competition Models Based on Block Growth Model [ J ]. Science and Technology Information, 2008 ( 34 ) : 4 + 6.

[2]Cotter, P. D., Ross, R. P., & Hill, C. (2012). Bacteriocins — a viable alternative to antibiotics? Nature Reviews Microbiology, 11(2), 95–105. doi: 10.1038/nrmicro2937