The modeling consists of three parts:
a) a modeling of the plane distribution of protons
b) the aggregation circuit of quorum sensing.

Part 1 The Plane Distribution of Proton

Theory and Method

we use the promoter PureI which is sensitive to the concentration of H+ to start the process of quorum sensing of E.coli. As H+ is the foundation of following steps, we need to build an environment with a proper gradient of protons' concentration. Therefore, we use a pole which we can put acid into, and take it as the maximum concentration of acid (Fig.1). We utilized two poles, one of which is instilled by acid that can supply H+. With the diffusion of H+, a gradient of protons' concentration is able to form. A schematic view of the specific equipment is revealed as following. Correspondingly, We build (stimulated) the spatial distribution of protons and get the solutions of diffusion according to Fick's Laws of Diffusion: $$C\left(r,t\right)=C_0\mathrm{erfc}\left(\frac{r}{2\sqrt{Dt}}\right)$$ $r$ is the radical distance from the pole; $t$ is the time;$C_0$ is the concentration of $H^+$ at the pole(that is the concentration of acid that we add);the length $2\sqrt{Dt}$ is the length of diffusion used to measure the length that $H^+$ diffuses in direction of x axis at time t. The length is called the diffusion length and provides a measure of how far the concentration has propagated in the x-direction by diffusion in time t .XX D is the diffusion coefficient; the diffusion coefficient of $H^+$ is $9.1\times10^{-8} m^2/s$, [1] http://www.ncbi.nlm.nih.gov/pubmed?term=15100022 erfc is complementary error function: $$\mathrm{erfc}=1-\mathrm{erf}(x)=\frac{2}{\sqrt{\pi}}\int_{x}^{\infty}e^{-t^2}\,dt$$

Simulation

Take the first two terms of the complementary error function after Taylor series expansion as the approximate number:
As a quick approximation of the error function, the first 2 terms of the Taylor series can be used: $$C\left(r,t\right)=C_0\left(1-\frac{x}{\sqrt{Dt\pi}}\right)$$ We can get the image as followed: fig.2 The plane distribution of protons in the culture medium
We divide the wasted water into two parts, the acid one and the alkaline one. Since the diffusion channel is quite small compared to the source of wasted water, we will not take the channel into our consideration. Thus, under the assumption of pH=6 in the acid part and pH=8 in the other part, and taking the source into consideration, we can get a new equation when t tends to infinitely great: $$C(r,\infty)=\lim_{t \to \infty}C_0\left(1-\frac{r}{\sqrt{DT\pi}}\right)=C_0$$ Therefore, when equilibrium is reached, pH in the acid part is always 6 and pH in the other part is 8. When these two parts neturalize after diffusion, the concentration of H+ decreases rapidly. image Fig. 3 the distribution of H+ in the medium

Part 2 The Aggregation Circuit of Quorum-sensing

Basic Concepts and the Aggregation Circuit

The ODE formalism models the concentrations of RNAs, proteins, and other molecules by time-dependent variables with values contained in the set of nonnegative real numbers. Regulatory interactions take the form of functional and differential relations between the concentration variables. For a typical transcription-translation process, the ODEs modeling approach associates two ODEs with any given gene $i$; one modeling the rate of change of the concentration of the transcribed mRNA $r_i$, and the other describing the rate of change of the concentration of its corresponding translated protein $p_i$. Thus for our network with 3 genes we have: $$\frac{dr_i}{dt}=F\left(R_i(p_1), R_i(p_2), R_i(p_3)\right)-\gamma_ir_i$$ $$\frac{dp_i}{dt}=P_i(r_i)-\delta_ip_i$$ Where (1) describes transcription, (2) describes translation, and $i = 1,…,N$. The functions $R_i(p_j)$ describe the dependence of mRNA concentration on protein concentration $p_j$ (If protein $p_j$ has no effect on mRNA $r_i$, then correspond function is set to zero.) The functional $F(•$) in (1) is defined in terms of sums and products of functions $R_i$. Function $P_i$ in (2) describes the translation of the mRNA $r_i$ into a protein $p_i$. Parameters $\gamma_i$, $\delta_i (i = 1,…,N)$, represent the degradation parameters of the mRNAs and proteins produced by gene i. As is common, we shall assume that the degradation of proteins or mRNAs is not regulated, namely that it does not depend on the concentrations of other molecules in the cell. Function $R_i$ is assumed to be in the form of Hill function as usual (since our cases are all inhibitors, we shall denote the Hill function $h^-(p,K,n))$, and the function $P_i$ is taken to be a linear term proportional to the concentration of mRNA $r_i$. $$h^-\left(p_i;K_i,n_i\right)=\frac{K_i^{n_i}}{p_i^{n_i}+K_i^{n_i}}$$ where $K_i$ is the microscopic dissociation constant, and $n_i$ is Hill coefficient, describing cooperativity. The regulatory network:

Equations and Parameters

Based on the fundamental assumptions above, we use ODEs to stimulate the circuit above. Taking both time and space into account, we can get all the equations as followed: (r represents the radical distance of one point in this environment from the pole.) \begin{equation} \frac{d[H^+]}{dt}=\frac{\delta C(r,t)}{\delta t} \end{equation} \begin{equation} \frac{d[r_1]}{dt}=\frac{m_1}{1+\left(\frac{[H^+]}{K_1}\right)^n}-\gamma[r_1] \end{equation} \begin{equation} \frac{d[LuxI]}{dt}=k[r_1]\delta_1[LuxI] \end{equation} \begin{equation} \frac{d[AHL]}{dt}=\frac{V_1[LuxI]}{K_2+[LuxI]}-k_1[LuxR][AHL]-\delta _2[AHL] \end{equation} \begin{equation} \frac{d[LuxR]}{dt}=k[r_1]-k_1[LuxR][AHL]-\delta_3[LuxR] \end{equation} \begin{equation} \frac{d[A]}{dt}=k_1[LuxR][AHL]-\frac{m_2}{1+\left(\frac{[A]}{K_3}\right)^n}-\delta_4[A] \end{equation} \begin{equation} \frac{d[cI]}{dt}=k[r_2]+\alpha-\delta_5[cI] \end{equation} \begin{equation} \frac{d[r_3]}{dt}=\frac{m_3}{1+\left(\frac{[cI]}{K_4}\right)}-\gamma[r_3] \end{equation} \begin{equation} \frac{d[LacI]}{dt}=k[r_3]=k[IPTG]-\delta_6[LacI] \end{equation} \begin{equation} \frac{d[r_4]}{dt}=\frac{m_4}{1+\left(\frac{K_5}{[LacI]}\right)^n}-\gamma[r_4] \end{equation} \begin{equation} \frac{d[CsgD]}{dt}=k[r_4]-\delta_7[CsgD] \end{equation} \begin{equation} \frac{d[CsgD]}{dt}=k[r_4]-\delta_7[CsgD] \end{equation} The parameter in these equations are listed in the table 1: Table 1

Meaning	Name	Unit	Value

For a certain depth of the biofilm, the concentration of oxygen is a constant in our model . Therefore we could solve these equations at different oxygen concentration and combine all the results to show how this system work. Because of the nonlinearity of the Hill functions, the solutions of a system of ordinary differential equations of a network of many genes cannot generally be determined by analytical means. http://www.ncbi.nlm.nih.gov/pubmed?term=15100022