Our modeling work is aimed at verifying the feasibility of our ternary symbiotic fermentation system, putting forward suggestions to the design of system by iteration, and finally simulating the expected behaviors.

On the one hand, we established mathematical models about the three microorganism mainly based on ODEs. First, we modeled the growth of Microbial Population and fitted our experimental data, the parameter values from which are used in future models. Then we simulated some gene pathways in E.coli to confirm the positive role of SacC in accelerating the sucrose utilization in E.coli. We also provided assistance to starvation promoter selection by simulating the pathways of different promoters. With the preliminaries above, we gave an insight into the quorum sensing pathways. The initial results are not ideal, hence we proposed a suitable degradation rate of AHLs and the necessity of improving sucrose production. We compared two schemes for regulating sucrose production and chose the better one.

On the other hand, we established a physical model about our specially designed fixed bed fermenter. It rendered a service to hardware design, including the location of entrances and exits, proportion of bacteria, arrangement of fixed beds and the shape of fermenter. Subsequently, we characterized the physical properties inside the fermenter, in order to verify its functions. Finally, we predicted the structure of optically controlled adhesion protein to calculated the adhesive force, which is compared with the fluid flush on bacteria.

All our modeling jobs are based on reasonable assumptions and reliable parameters, interacting well with wet lab.

We have built mathematical models to predict the growth of three kinds of microorganism in our system. In order to capture the effect of limiting substrate concentration, we chose the Monod Equation

which was put forward by Monod in 1942 and has the same form as Michaelis-Menten Equation. The parameter `\mu_{max}` represents the maximum growth rate, `[S]` represents limiting substrate concentration, and `K_{S}` is the substrate concentration at which the growth rate reaches 1/2`\mu_{max}`.

To visualize the modeling, we tentatively assume that `\mu_{max}` is 0.05 min^-1, `K_{S}` is 9e8 `\mu g\cdot L^{-1}`, and `[S]` is 1.5e6 `\mu M`. The initial quantity of bacteria is set as 1e6, and S is regarded as glucose. Soluting this ODE by ode45 solver in MATLAB, the growth curve for bacteria given by this model is obtained as below.

However, this model hasn't consider the influence of environmental capacity, namely the maximum population quantity. So we revised the model by combining Logistic equation

with Monod equation. The hybrid Logistic-Monod equation

is obtained by multiplying the Monod equation with the self-inhibiting factor `(1-N/N_{m})`.

We soluted the new ODE when `N_{m}` takes different values such as 1e8, 5e8, 1e9 and 2e9. The result shows that the population quantity will come to a stationary phase in the end, and the final quantity is determined by `N_{m}`. Then we fixed `N_{m}` as 2.1e8, and plotted a three-dimensional image to show the variation of growth curve depending on the concentration of limiting substrate.

However, there are still defects in the model, for we hasn't consider the decline of bacteria population caused by depletion of nutrients and accumulation of harmful substances. Therefore, we ulteriorly revised the model by introducing a coefficient `d` which represents the attenuation rate of population quantity. Thus the revised hybrid Logistic-Monod equation is

We soluted the final ODE when `N_{m}` is 2.1e8 and `d` takes different values such as 1e-4, 2e-4, 3e-4 and 4e-4. The result shows the population quantity will decline after a certain length of stationary phase, and the attenuation comes early as `d` rises. Then, as above, we fixed `N_{m}` as 2.1e8, and plotted a three-dimensional image to show the variation of growth curve depending on the concentration of limiting substrate.

Finally, we can give the Revised-Double-Substrates-Logistic-Monod equation

The equation can be written as the ODE above only when `[NH_{4}^{+}]` and `[Glucose]` are considered as constant. Under the circumstances, we can give the analytical solution for this ODE:

Next step, we fitted the growth curve obtained in wet lab experiments using the revised Logistic-Monod equation above, from which we can acquire the parameters being used in the future modeling. The concentration of glucose was assumed as 1.5e6 `\mu M`, which is considered appropriate for E.coli in reference[1]. Below are the fitted curves and values of parameters for E.coli and A.caulinodans.

Being a kind of photoautotroph, S.elongatus will grow under the limiting condition of light intensity. Thus we can regard the light intensity `I`(or number of photons) as the limiting substrate `[S]` in Monod equation

Considering that light is continuously absorbed as passing through the culture solution, the intensity varies with the position in the conical flask. This is called "self-shading"[2].

According to Beer-Lambert Law

where A is absorbance, T is transmittance, `I_{0}` is incident intensity, `I` is emergent intensity, `b` is thickness of Absorption layer, `c` is molarity of light-absorbing substance and `K_{m}` is molar absorption coeffcient. Because we didn't find exact values for all these parameters, this model is only a qualitative one.

We regard the conical flask as a cylinder with a radius of 3.5 centimeters, and calculate the intensity distribution. The cylinder is similar to the shape of our fermenter, which will be well simulated. Here we used absorption coefficient instead of molar absorption coefficient

Plot this function:

Plug `I(r)` into Monod equation, then the growth model of the entirety is expressed as an integral over the position

This resembles setting weights for cells exposing to light of different intensities, which will make the S-shaped curve tend to be a straight line. And the maximum population quantity will also decline. The growth curve before and after considering light distribution are shown below.

This model gives an insight into the growth of S.elongatus on the fixed bed inside the fermenter. Afterwards, we excitedly found that our experimental results are consistent with this model.

In this model, we aim to verify the effectiveness of the plasmid pUC-SacC and quantitatively calculate how SacC helps E.coli intake sucrose as a validation model for the wet lab.

We only focus on the changes of maximum growth rate r in Logistic equation influenced by sucrose, so we use the hybrid Logistic-Monod equation above to research.

According to the result of the experiment part, we can plot the growth curves of E.coli in the sucrose medium with respectively pUC-empty and pUC-SacC:

On the basis of the logistic equation, we can solve out the general model to describe the growth curve of E.coli as:

By setting the number of OD600 in the plateau phase as the maximum growth limit of E.coli, we may fit the curve of the model with the given data.

The growth curves are

Therefore, we can get that the growth rate of E.coli with these two different plasmids are separately 0.5504 (`h^{-1}`) and 0.1607(`h^{-1}`).

By searching and reading papers, we found a study showing that the maximum growth rate of E.coli under 37°C without any limited substrates is approximately 0.54[3], which is rather close to the calculated growth rate in our model. Therefore, we can assert that the presence of SacC makes sucrose no longer a limited substrate in the medium, verifying the effectiveness of SacC.

Besides, we introduce the Monod equation to describe the growth rate of E.coli under limited substrates:

In which `\mu_{max}` is 0.5504(`h^{-1}`), the growth rate of E.coli with pUC-SacC. And from the reference, we got to know that the performance of E.coli under a high concentration of sucrose is like that of under 0.2% of sucrose[4]. So according to the data in the paper, we can set the S in our model as 2g/L to calculate Ks which is 4.85 g/L. By comparing the Ks in the model with that in the reference[5], we can know that the Ks in the model is reasonable. Therefore, we can verify the correctness of our model. Moreover, we can quantitively work out that the presence of SacC can increase the growth rate of E.coli by 3.42 times to an unlimited rate against sucrose.

This model acts as a validation model, confirming the correctness of the data from the wet lab, verifying the effectiveness of pUC-SacC in improving the efficiency of E.coli to intake sucrose, and quantitatively determining the effect of SacC as 3.42 times’ improvement.

This model set up before the wet lab concerning carbon starvation promoter is made to pre-confirm the assumption that lowering the carbon concentration in the medium can increase the promoter strength of carbon starvation promoter in E.coli.

By searching references, we noticed a substance called cAMP-CRP complex, which was found directly correlated to the increase of promoter strength under starvation [6][7]. The reference[6] told us that when the starvation signal is input, the activation network of the E.coli will synthesize the cAMP-CRP complex, and reference[7] told us the concentration of cAMP-CRP has a strong positive correlation with the strength of starvation promoter in a certain range of concentration.

Therefore, by simulating and solving the activation network of cAMP-CRP, we can forecast the increase of promoter strength under starvation.

According to the activation network, we can list the ODEs for it:

Due to the fact that hardly can all of the parameters in the ODEs be clearly found, and for many parameters, only their order of magnitude is ascertained. So after estimating、 concerning the simple relationship between parameters (for example the difference between orders of magnitude is often 1 or 2 between k of the forward and backward reaction). So we determine the parameters `k_{1}`=1e3, `k_{-1}`=1e2, `k_{2}`=3.3e2, `k_{4}`=1e2, `k_{-4}`=1e1, which are exact, and `K_{1}`, `K_{2}`, `K_{3}`, `K_{4}`, are the general order of magnitude of rate from which the gene generate protein, -1`\sim`2, and therefore, we can solve the ODEs with the initial value of carbon source at 0.8M.

Therefore, we can justify that the carbon starvation do strengthen the promoter. So, the wet lab’s premise that carbon starvation can improve the strength of promoter is verified.

The quorum sensing system we designed has two orthogonal feedback circuits. The starvation promoters are activated to express CinI or RhlI when E.coli is deficient in carbon or nitrogen sources, whereupon the signal moleculars AHL1 or AHL2 are produced and diffuse to S.elongatus and A.caulinodans. In the latter two bacteria, there are constitutive promoters expressing the signal receivers CinR or RhlR, which can form heterotetramers with AHLs to activate genes involved in nutrient production. After nutrients are accumulated sufficiently, AHLs production stops, preventing overproduction of nutrients. Consequently, it is expected that the concentrations of nutrients and AHLs form an oscillation while the E.coli population maintains at a certain level.

Considering the two quorum sensing pathways are orthogonal, we simulated them respectively to avoid a complicated ODE system. The reactions and ODEs in Three kinds of bacteria are introduced separately below.

Signal moleculars are named as AHLi for convenience. The corresponding synthetases are Syni and receivers are Ri (where i represents 1, 2). Then the reactions related to the expression of Syni are

The rates of the latter three reactions follow the Law of Mass Action, where `k_{tl}` is translation rate, `d_{mRNA}` is degradation rate of mRNA, and `d_{pr}` is degradation rate of protein. The transcription rate `k_{tr}` is described by Hill function `k_{tr}=a_{0}+a\cdot \frac{K^{n} }{K^{n}+[S]^{n}}`, where `a_{0}` is basal transcription rate, `a` is maximum transcription rate, `K` is repression coefficient, `[S]` represents the concentration of carbon or nitrogen source. Thus we can give the related ODEs

To simplify the modeling, we approximated Hill function using logic input function, namely a step-function `k_{tr}=a_{0}+a\cdot \text{heaviside}([S_{critical}]-[S])`.

The biochemical reaction circuit of AHL synthesis is drawn as below,

which can be simplified as

This is a double-substrate enzymatic reaction following the ping-pong mechanism. So the reaction rate can be given by Double-Substrate Michaelis-Menten Equation

where `K_{cat}` is catalytic number and `K_{M,i}` is Michaelis constant of `S_{i}` (i=1,2). Now we can give the ODEs

where `k_{SAM}, k_{aACP}, d_{SAM}, d_{aACP}` are respectively produce and decay rate of SAM and aACP. `AHL_{E}` represents the AHL concentration in E.coli while `AHL_{e}` represents that in environment. The transmembrane diffussion rate `D` is calculated by `D=\frac{S\cdot P}{V}` where `S` is the surface area of bacterium, `P` is the membrane permeability, `V` is the volume of bacterium.

The growth function of E.coli

is given by the previous model, the parameters of which are listed in Table.1.

Several proteins, including `R1`, `\text{nifA}`, `RpoN(\sigma_{54})`, `NIF`, are expressed in A.caulinodans. Their reactions are similar to the ones of Syni expression shown above. But there are some differences among their reaction rates. Firsty, we can write the ODEs about R2 expression,

where `k_{tr,c}` represents the transcription rate of constitutive promoter.

Before calculating other rates, let's model the process of RpoN activating the expression of NIF. As we all know, there are some kinds of `\sigma`-factors in bacteria can be associated with core enzymes(represented by "RNAP") to form holoenzymes. Different `\sigma`-factors compete for the limited core enzymes to activate different promoters. We divided all kinds of `\sigma`-factors into `\sigma_{54}`(RpoN) and `\sigma_{i}`(any other `\sigma`-factors). Their association and dissociation rates with RNAP `k_{54}`, `k_{-54}`, `k_{i}`, `k_{-i}` are consulted in reference[12]. Now we can give the reactions of holoenzymes production:

and the reaction of transcription with RNAP.`\sigma_{54}` taken into account:

Now we can have the ODEs about coreenzymes and holoenzymes

and the ones about NIF expression

Holoenzyme quantity influences `k_{tr}` linearly, while the effect of nifA follows Hill function.

When reaching A.caulinodans, AHLs produced by E.coli will go through another transmembrane diffusion. Then they associate with R1 followed by a dimerization. Related reactions and ODEs are as follows:

where `AHL_{A}` represents AHLs in A.caulinodans, `AHL_{e}` represents AHLs in environment, and `k_{1}`, `k_{-1}`, `k_{di,1}`, `k_{di,-1}` are association and dissociation rates.

Finally, the heterotetramer `(AHL1_{A}:R1)_{2}` can activate the expression of nifA and RpoN, which gives the ODEs

Related parameters are acquired by fitting experimental datas in wet lab. Now, we can write the ODEs about `\sigma`-factors.

The growth function of A.caulinodans

is given by the previous model, the parameters of which are listed in Table.2.

Finally, we simplified the regulation of nitrogen production. It was considered that NIF could promote producing ammonium and the specific procedures are ignored. The promotion roughly follows Hill function. Suppose that the basal ammonium produce rate is `k_{n}`, the maximum rate after activated is `k_{nm}`, and the activation coefficient is `K_{NIF}`. Because of the influences of other limiting conditions in A.caulinodans, the acceleration of ammonium production will be remarkable but come to a threshold later. Therefore, we think the hypothesis to be reasonable. Then we introduced the consumption rate of ammonium `d_{NH_{4}^{+}}`. Both the production and consumption are related to the quantity of cells, so we multiplied them with the factors `\frac{N_{2}V_{cell}}{V_{e}}` and `\frac{(N_{1}+N_{2})V_{cell}}{V_{e}}`. Now we have the ODE of ammonium:

Before soluting the ODE set above, we inspect it and found that there is no degradation mechanism for AHL, which may lead to malfunction of feedback circuits. This is also mentioned in reference. We soluted the system of ODEs using the ode15s solver in MATLAB, and confirmed this inference.

According to Fig.18 , ammonium is consumed untill less than critical value. Then it will be called back by AHL. But it won't decline again, for AHL can't be degradated. This will cause the surplus of nitrogen source, which is not energy friendly. Therefore, we appended degradation rates of AHLs to equations of AHL concentration. The degradation mechanism can be realized by introducing AHL degrading enzymes or adding degradation tag to Syni.

Solving the ODE set by ode15s again, and iterating the degradation rate of AHL, we obtained an ideal result.

As shown in the figures above, ammonium is heavily consumed during logarithmic phase of E.coli frowth. When ammonium decreases to the critical value, AHL will proliferate to accelerate ammonium production. And this feedback process will repeat over and over, forming a steady oscillation. With ammonium concentration maintaining close to critical value, E.coli can keep a certain population quantity. The estimated degradation rate of AHL is around 0.1`\mu M\cdot \text{min}^{-1}`, which gives guidance to wet lab experiments. Besides, the maximum transcription rate of starvation promoter is required to be around 200 times higher than the initial rate, which needs to be implemented in experiments.

Proteins expressed in S.elongatus are `R2` and `cscB`. The related reactions are similar to the ones of Syni. And the reactions about AHL2 are also similar to the ones about AHL1. Thus we can write the ODE set below, where we had added into the degradation rate of AHL.

In the above equations, we simplified the processes of sucrose production in S.elongatus and degradation in E.coli. Instead, we regard hexoses as the carbon source, and cscB plays a role in accelerating its transmembrane transportation. However, the solution shows that the increase of hexose production can't afford the growth of E.coli. Concentration of hexose can't go back upon the critical value, so E.coli will starve to death in the end.

Therefore, it is demanded to improve our system design. We suggested appending a mechanism for increasing sucrose production in S.elongatus in addition to increasing cscB. Then we iterated the parameter `k_{nm}` and solved the ODE set again. The results show that when sucrose(or hexose) production is 2~3 times higher, there will be oscillation behaviors similar to that in the previous model.

In conclusion, we simulated and iterated the quorum sensing circuits. Some suggestions for improvement were given. And the improved model verified the feasibility of our system, namely the three microorganism can form a stable symbiotic relationship to achieve our production objectives.

All the related parameters are listed in the end.

In the Quorum Sensing model, we found that regulating cscB cannot fit our demand of the concentration of sucrose. So we need to decide new ways for S.elongatus to produce more sucrose.

In order to predict the rise of sucrose production, we built the carbon assimilation pathway kinetic model of S.elongatus and simulated this process. Below is the biochemical reaction diagram of S.elongatus.

Because we only focus on the changes of relative concentration between sucrose and glycogen, the concentration of F6P and the former reactants are not important. Therefore we might as well suppose F6P has a fixed concentration. Then the critical reactions of this pathway are as follows:

ODEs:

We build the model in Matlab Simbiology. The diagram of this model is shown here.

We anticipate two ways of gene silencing to improve the production of sucrose. One is to directly silence the gene of producing glycogen (the glgA gene) to decrease its direct competition against sucrose production. The other is to silence the gene of producing sucrose (the spp gene) to accumulate glycogen, later when the signal of starvation occurs, we remove the silencing of spp and silence glgA alternatively. The cumulated glycogen will improve the generation rate of sucrose to a high degree, therefore making the S.elongatus produce more sucrose in the same time.

We first tried silencing glgA and get the following image. (the different colours of curves are different groups of initial value of AGPase and SSase).

We can see that after gene silencing of glycogen gene glgA, the production of sucrose became about 3.3e4 micromolarity at 1800 minute.

Then we try to silence sucrose gene spp within the first 900 minutes, and then remove the gene silencing and get the following image:

We can see that at 1800 minute, the production of sucrose became about 3.7e4 micromolarity, which is about 10% improvement compared with former method.

This model proves that firstly silence spp gene and then remove that can improve the production of sucrose to some degree. Even though the increase of production in reality may not be so big due to the limit of enzyme activity maximum or the limit of transfer process, this model can still guide the wet lab to adopt a new approach to improve the production of sucrose.

Unlike traditional fermenters, each of our fermenters is connected to the other two, and there is interaction between them, which means that the flow velocity distribution in the fermenter should be as uniform as possible to ensure that there is no retention of nutrients and signaling molecules.

Here we tested each case using Ansys Fluent. Fluent solves fluid problems by finite-volume method, with internal nodes of an unstructured grid placed in the flow field in an irregular manner (more elaborate in complex structures and looser in regular areas). By transforming the equations of continuity equation, momentum, energy, and so on, into universal forms, we got:

In this equation, `\phi` is generalized variables, which can be velocity, temperature or concentration of some of the physical quantities to be sought. `\Gamma` is corresponding generalized diffusion coefficient of `\phi`, `S` is the generalized source term. The boundary value of the variable `\phi` in the control body is known. Then we integrated the volume of the control body at the node, and we got:

`\Delta V` is the volume value of the control volume. When the control volume is very small, `\Delta V` can be expressed as `\Delta V\cdot A`, where A is the area of the control volume interface. Thus:

If the mesh is uniform, do linear interpolation:

Then we will have:

And the we let `S=S_{C}+S_{P}\phi_{P}` (for the source term S, it is usually a function of time and physical quantities) , then at each node there is a discrete form:

Fluent can get the value of each physical quantity at each node by solving the equations.

Note bene: Unless otherwise noted, the algorithm used in Fluent by default is the Coupling method (there are also methods like SIMPLE, SIMPLEC, PISO, and so on in Fluent).

The coupling method is used to solve discrete equations at the same time, and the variables (et al.) are solved simultaneously. The solving process is as follows.

(1) the initial pressure and velocity are assumed to determine the coefficients and constants of the discrete equations.

(2) solving continuous equation, momentum equation and energy equation simultaneously.

(3) solving turbulence equations and other scalar equations.

(4) to judge whether the calculation on the current time step is convergent or not. If convergent, return to the second step, iterative calculation; if convergent, repeat the above steps, and calculate the next time step of the physical quantity.

For the traditional fermenter, which is mainly composed of porous media and fixed phase, the flow velocity distribution affected by the entrance location is only in the area near the entrance because of the huge resistance inside. Figures are shown below:

The flow velocity at the entrance of Figure 31 and 32 is 1 m/s, and that at the entrance of Figure 33 is 0.5 m/s. The whole region is set as porous medium, fluid as water. The gravity downward is 4.9 `m\text{/}s^{2}` (Mars gravity). The other parameters are all default.

Figure 31, Figure 32 and Figure 33 respectively shows the flow velocity distribution at different inlet and outlet positions with the flow rate unchanged. It can be seen that the flow in the porous medium is mainly influenced by gravity when the inlet velocity is small, so it is quite uniform away from the inlet and outlet. At the entrance and exit, the flow velocity can be reduced by increasing the cross-section, that is to say, the area of non-uniform flow velocity can be reduced without changing the flow rate, this is actually consistent with our experience with fermenters (an entrance, but there is a gap below the entrance where there is no filler-equivalent to an entrance the size of the tank) .

However, the design of the fermenter for S.elongatus in our project is obviously not suitable to maintain this level of density due to the need for light, additional simulations were performed for fermenters with sparse stationary phases (as opposed to infiltration, where different exit locations would result in a larger flow rate distribution fluctuation) :

The velocity at the entrance is 1m/s except for 0.5 m/s in Figure 32. The whole region fluid is water, gravity down 4.9 m/s 2, the other parameters are all default. As you can see from Figures 31 and 32, while keeping the traffic constant, for such a sparse fermenter, the overall flow rate can not be made more uniform by increasing the inlet area and decreasing the inlet flow rate; As can be seen from Figure 34 and 35, in the case of sparse stationary phase, the change of one inlet can cause a large change of velocity distribution.

Figure 35, Figure 36 and Figure 37 are better than the others on the condition that the span and the minimum velocity distribution are as large as possible, and Figure 35 is the optimal distribution of inlet and outlet.

Our project involves the collaboration of multiple bacteria, some of which act as carbon sources for symbiotic systems and others as nitrogen sources, while each type of bacteria continues to consume nutrients. We want to maximize our productivity by maximizing the number of bacteria that can produce the products we need so that the system doesn't have a nutrient surplus.

According to the production rate and consumption rate of each bacterium, the number of S.elongatus was assumed to be x, the number of A.caulinodans to be y, and the number of E.coli to be z. Then we can list the following mathematical formula (inequality 3 is set up because we are only trying to find the ratio of bacteria to bacteria in a given population) :

This kind of problem can be solved simply by one-way method (when the constraint condition is only inequality condition and the objective function is linear, the optimal proportion can be solved by one-way method). The corresponding ratio of x, y and z is 125:120:55.

In order to avoid the influence of transportation and uneven connections between the fermenters, and to keep the expected ratio. The connections between the fermenters must be topologically circular and uniformly distributed, as shown in the following figure (the ratio here is approximately 2:2:1):

Of course, for conservative reasons, in order to increase the stability of the symbiotic system, the number of bacteria responsible for producing the nutrients should be raised to ensure that the system is not destabilized.

The fermenters simulated in the model above have existing fixed-bed structures, which may not satisfy the required efficiency on Mars. Therefore, further optimizations are in demand, for which we have come up with some ideas. We first selected the most suitable structure by self-iteration according to certain standards, and then used computational fluid dynamics simulation to verify whether the design meets our requirements.

We have summarized these jobs in a pdf file: Design of Fixed-bed in Fermenter

Bacteria living in a fermenter have a requirement for nutrient concentrations, and we need to ensure that at the right flow rate and entry concentration, the concentration of the bacteria on the surface of the whole fermenter is within the allowable range of their survival.

In order to simulate the behavior of secreting and consuming substances, we chose the component transport model in Fluent and the fermenter with dense heap (which the E. coli lives in) as the main simulation object , taking into account the previous simulation results (which showed that the flow velocity distribution was uniform and was most uniform when the inlet radius was the same as the radius of the fermenter) , we chose to arrange some tube walls in the fermenter as the outlet for releasing glucose at a fixed rate to simulate the state of glucose consumption by the bacteria. And based on previous calculations about the bacterial mix, considering that our system is not stable, the cyanobacteria mix will end up being higher, the amount of glucose entered per unit time is about twice the total glucose consumption in the fermenter, on this basis, the inlet and outlet flow velocity is calculated to be about 800 times of the precipitation flow velocity (at this time, the ratio of glucose in the precipitation flow velocity is 0.2 if it is too large, it means that the precipitation flow velocity is smaller when the total amount of the precipitation glucose remains constant, and if the grid segmentation is not precise enough, the precipitation flow velocity is too small and the concentration change caused by it will be ignored in the calculation. By adjusting the inlet flow rate and concentration to 800 times, the effect of simulated nutrient uptake on the overall flow rate distribution can be ignored (at the macro level) , and the effect on the concentration of the substance can not be ignored at the macro level) . The results were as follows:

It can be seen that the flow velocity at the entrance is 0.00316 m/s, the average internal velocity is 0.00312 m/s, and the inlet glucose concentration is 0.15 kg/L, the concentration in the fermenter met the glucose requirement range of 0.01-0.18 kg/l for bacteria.

According to the simulation results above, we confirmed the effectiveness of the fermenter system we designed.

During the project design, to achieve the replacement and plug-and-play nature of the strains in the fermenter, we designed a mechanism to produce colonies controlled by blue light. But time constraints only allowed this interesting design to be briefly verified by modeling.

In the mechanical modeling, we will predict the binding free energy between two mutated proteins. We use haddock2.4 for hohomology modeling and docking, get its predicted binding structure and binding free energy [13].

Considering that we wanted our bacteria to be artificially controlled to adhere to the fermenter, we selected a pair of light-controlled adhesion proteins called pMag and nMag. They were homologs of Vivid (VVD) from Neurospora crassa with a mutation at positions 52 and 55, respectively, and thus each carries a positive and negative charge [14].

Since there was no literature specifying the three-dimensional structure and binding free energy of the two proteins, we decided to predict the structure of the two proteins and dock them using homology modeling of swissmodel. Because the conformations of pMag and nMag possess the protein-protein interaction only after excitation by blue light irradiation, the template we selected was one of the VVD that was excited by blue light (PDB: 3d72) to become a dimer before the mutation [14].

By docking the pdb file obtained from the same source model using haddock2.4, and selecting positions 11-20 (i.e. alpha-helix at VVD47-56) as the docking region, we got seven predicted cluster.

By comparing the Z-score, Cluster Size and the electrostatic potential on the surface of docking of predicted structure, cluster 7 was the best whose

We selected a total of 20 amino acids on the alpha-helix of the two proteins at the binding site of this docking structure and treated the peptide bond as two positive and negative charges that converge by about 1Å. Based on the above, we calculated the interaction force between them, i.e., the orientation force, which was a repulsive force of about `-2.4 × 10^{-11}N` using the Coulomb force equation. We then calculated the forces between the four charged side chain amino acids to be about `-6.1 × 10^{-10}N`. The magnitude of the electrostatic interaction was about one order of magnitude larger than the orientation force, so it was known that the van der Waals force, which was in the same order of magnitude as the orientation force, does not produce large errors in the results.

Considering that the average value of flow velocity and the degree of fluctuation in the fermenter with dilute media was much larger than that in the fermenter with dense stack media, the flow velocity at the entrance of the dilute fermenter was set to 1m/s and the maximum flow velocity was about 2m/s, and the glucose concentration was about 15% (to ensure that the glucose mass fraction was close to 18%).

Based on this, we obtained the magnitude of the scouring force caused by the water flow on the bacteria. Based on this data, we performed the corresponding two-dimensional simulation of the spheres carrying bacteria in a fermenter using a dilute medium to obtain the maximum flow velocity gradient on the surface of the spheres. This value was used to calculate the maximum tangential stress, i.e. the maximum value of the scouring force of the fluid on the bacteria, using Newton's equation for viscous forces.

The maximum flow rate gradient was about `4×10^{-3}s^{-1}`, and the corresponding tangential pressure was `1.2×10^{-5}Pa` (considering that the viscosity of liquid glucose was about `3.7×10^{-3}Pa\cdot s` at room temperature, while that of water was `2.98×10^{-3}Pa\cdot s`, the difference was not large, so we directly take `3×10^{-3}Pa\cdot s` for estimation). ×Since it was densely distributed, we take half of its cross-sectional area as the corresponding area subject to scouring pressure. The corresponding scouring force was `4.6 × 10^{-16}N`, which was much smaller than the viscous force between proteins.

parameter	value	unit	description	reference
`a_{0}`	0.03	`\mu M\cdot \text{min}^{-1}`	basal transcription rate	assumption
`a_{0,NR}`	0.05	`\mu M\cdot \text{min}^{-1}`	basal transcription rate of nifA:RpoN	wetlab
`a_{0,NIF}`	0.0045	`\mu M\cdot \text{min}^{-1}`	basal transcription rate of nifA:RpoN	assumption
`a_{NR}`	27	`\mu M\cdot \text{min}^{-1}`	maximun transcription rate of nifA:RpoN	wetlab
`a`	29.97	`\mu M\cdot \text{min}^{-1}`	maximum transcription rate	assumption
`K_{M,NR}`	1.67e-5	`\mu M`	activation coefficient of nifA:RpoN transcription	wetlab
`a_{Syn1}`	20	`\mu M\cdot \text{min}^{-1}`	maximum transcription rate of Syn1	assumption
`n_{NR}`	1.107	1	Hill coefficient of nifA:RpoN transcription	wetlab
`k_{tl}`	10	`\text{min}^{-1}`	translation rate	assumption
`k_{tl,Syn1}`	20	`\text{min}^{-1}`	translation rate of Syn1	assumption
`k_{tl,NR}`	9	`\text{min}^{-1}`	translation rate of nifA:RpoN	assumption
`d_{mRNA}`	0.03	`\text{min}^{-1}`	degradation rate of mRNAs	[8]
`d_{pr}`	0.0033	`\text{min}^{-1}`	degradation rate of ptoteins	assumption
`d_{pr,NR}`	0.33	`\text{min}^{-1}`	degradation rate of nifA:RpoN mRNA	assumption
`k_{SAM}`	51	`\mu M\cdot \text{min}^{-1}`	produce rate of SAM	[12]
`d_{SAM}`	0.01	`\text{min}^{-1}`	degradation rate of SAM	[12]
`k_{aACP}`	50	`\mu M\cdot \text{min}^{-1}`	produce rate of aACP	[12]
`d_{aACP}`	0.01	`\text{min}^{-1}`	degradation rate of aACP	[12]
`K_{cat,1}`	0.46	`s^{-1}`	catalytic number of Syn1	[9]
`K_{M,SAM}`	14	`\mu M`	Michaelis constant of Syn1 for SAM	[10]
`K_{M,aACP}`	6	`\mu M`	Michaelis constant of Syn1 for aACP	[10]
`D`	0.033	`\text{min}^{-1}`	diffusion coefficient of AHL	wiki of Ecuador 2021
`V_{cell}`	1e-15	`L`	volume of bacterial cell	assumption
`V_{e}`	2e-3	`L`	volume of environment	wiki of Ecuador 2021
`\mu_{max,1} `	0.05019	`\text{min}^{-1}`	maximun growth rate of E.coli	wetlab
`\mu_{max,2} `	0.09	`\text{min}^{-1}`	maximun growth rate of A.caulinodans	wetlab
`k_{tr,c}`	0.1	`\text{min}^{-1}`	transcription rate of constitutive promoter	assumption
`k_{n}`	2.5e3	`\text{min}^{-1}`	basal produce rate of ammonium	assumption
`k_{nm}`	2.5e3	`\text{min}^{-1}`	maximum produce rate of ammonium	assumption
`k_{1}`	0.35	`\mu M^{-1}\cdot \text{min}^{-1}`	association rate of AHL1 and R1	[12]
`k_{-1}`	2	`\text{min}^{-1}`	dissociation of AHL1 and R1	[12]
`k_{di,1}`	0.03	`\mu M^{-1}\cdot \text{min}^{-1}`	dimerization rate of AHL1:R1	[12]
`k_{di,-1}`	3	`\text{min}^{-1}`	depolymerization rate of `(AHL2:R2)_{2}`	[12]
`N_{m,1}`	2.1e8	cell	maximum population quantity of E.coli	wetlab
`N_{m,2}`	8.5e7	cell	maximum population quantity of A.caulinodans	wetlab
`N_{0,1}`	1e6	cell	initial population quantity of E.coli	wetlab
`N_{0,2}`	5e5	cell	initial population quantity of A.caulinodans	wetlab
`K_{N}`	5.16e6	`\mu g\cdot L^{-1}`	activation coefficient of growth by ammonium	wetlab
`d_{1}`	1e-6	`\text{min}^{-1}`	decay rate of E.coli population	wetlab
`d_{2}`	8.2e-6	`\text{min}^{-1}`	decay rate of A.caulinodans population	wetlab
`k_{54}`	1	`\mu M^{-1}\cdot \text{min}^{-1}`	association rate of `\sigma_{54}` and RNAP	[11]
`k_{-54}`	0.3	`\text{min}^{-1}`	dissociation rate of `RNAP.\sigma_{54}`	[11]
`k_{i}`	1	`\mu M^{-1}\cdot \text{min}^{-1}`	association rate of `\sigma_{i}` and RNAP	[11]
`k_{-i}`	0.33	`\mu M^{-1}\cdot \text{min}^{-1}`	dissociation rate of `RNAP.\sigma_{i}`	[11]
`k`	1e8	`\text{min}^{-1}`	rate constant of NIF transcription	assumption
`K_{nifA}`	500	`\mu M`	activation coefficient of NIF transcription	assumption
`K_{NIF}`	500	`\mu M`	activation coefficient of ammonium produce	assumption
`d_{AHL1}`	0.1	`\mu M\cdot \text{min}^{-1}`	degradation rate of AHL1	assumption
`d_{NH_{4}^{+}}`	2.05e5	`\mu M\cdot \text{min}^{-1}`	consumption rate of ammonium	assumption

parameter	value	unit	description	reference
`a_{0,Syn2}`	0.03	`\mu M\cdot \text{min}^{-1}`	basal transcription rate of Syn2	assumption
`a_{Syn2}`	20	`\mu M\cdot \text{min}^{-1}`	maximum transcription rate of Syn2	assumption
`a_{0,sRNA}`	0.8	`\mu M\cdot \text{min}^{-1}`	basal transcription rate of sRNA	assumption
`a_{0,Hfq}`	0.8	`\mu M\cdot \text{min}^{-1}`	basal transcription rate of Hfq	assumption
`a_{sRNA}`	25	`\mu M\cdot \text{min}^{-1}`	maximum transcription rate of sRNA	assumption
`a_{Hfq}`	25	`\mu M\cdot \text{min}^{-1}`	maximum transcription rate of Hfq	assumption
`a`	29.97	`\mu M\cdot \text{min}^{-1}`	maximum transcription rate	assumption
`K_{M}`	40	`\mu M`	activation coefficient	assumption
`n`	2	1	Hill coefficient	assumption
`k_{tl}`	10	`\text{min}^{-1}`	translation rate	assumption
`d_{mRNA}`	0.03	`\text{min}^{-1}`	degradation rate of mRNAs	[8]
`d_{pr}`	0.0033	`\text{min}^{-1}`	degradation rate of proteins	assumption
`d_{Hfq}`	0.033	`\text{min}^{-1}`	degradation rate of Hfq	assumption
`k_{SAM}`	51	`\mu M\cdot \text{min}^{-1}`	produce rate of SAM	[12]
`d_{SAM}`	0.01	`\text{min}^{-1}`	degradation rate of SAM	[12]
`k_{aACP}`	50	`\mu M\cdot \text{min}^{-1}`	produce rate of aACP	[12]
`d_{aACP}`	0.01	`\text{min}^{-1}`	degradation rate of aACP	[12]
`K_{cat,2}`	0.46	`s^{-1}`	catalytic number of Syn2	[9]
`K_{M,SAM}`	14	`\mu M`	Michaelis constant of Syn2 for SAM	[10]
`K_{M,aACP}`	6	`\mu M`	Michaelis constant of Syn2 for aACP	[10]
`D`	0.033	`\text{min}^{-1}`	diffusion coefficient of AHL	wiki of Ecuador 2021
`V_{cell}`	1e-15	`L`	volume of bacterial cell	assumption
`V_{e}`	2e-3	`L`	volume of environment	wiki of Ecuador 2021
`\mu_{max,1} `	0.05019	`\text{min}^{-1}`	maximum growth rate of E.coli	wetlab
`\mu_{max,3} `	0.0004	`\text{min}^{-1}`	maximun growth rate of S.elongatus	wetlab
`k_{tr,c}`	0.1	`\text{min}^{-1}`	transcription rate of constitutive promoter	assumption
`k_{c}`	2.5e4	`\mu M\cdot \text{min}^{-1}`	basal produce rate of glucose	assumption
`k_{cm}`	2.5e4	`\mu M\cdot \text{min}^{-1}`	maximum produce rate of glucose	assumption
`k_{2}`	0.35	`\mu M^{-1}\cdot \text{min}^{-1}`	association rate of AHL2 and R2	[12]
`k_{-2}`	2	`\text{min}^{-1}`	dissociation of AHL2 and R2	[12]
`k_{di,2}`	0.03	`\mu M^{-1}\cdot \text{min}^{-1}`	dimerization rate of AHL2:R2	[12]
`k_{di,-2}`	3	`\text{min}^{-1}`	depolymerization rate of `(AHL2:R2)_{2}`	[12]
`N_{m,1}`	2.1e8	cell	maximum population quantity of E.coli	wetlab
`N_{m,3}`	2.1e8	cell	maximum population quantity of S.elongatus	wetlab
`N_{0,1}`	1e6	cell	initial population quantity of E.coli	wetlab
`N_{0,3}`	1e6	cell	initial population quantity of S.elongatus	wetlab
`K_{C}`	9e8	`\mu g\cdot L^{-1}`	activation coefficient of growth by glucose	wetlab
`d_{1}`	2e-6	`\text{min}^{-1}`	decay rate of E.coli population	wetlab
`d_{3}`	2e-6	`\text{min}^{-1}`	decay rate of S.elongatus population	wetlab
`K_{Hfq}`	350	`\mu M`	activation coefficient of glucose produce	assumption
`d_{glc}`	2e6	`\mu M\cdot \text{min}^{-1}`	consumption rate of glucose	assumption
`d_{AHL2}`	0.1	`\mu M\cdot \text{min}^{-1}`	degradation rate of AHL2	assumption

parameter	value	unit
`V_{max,GPI,f}`	317200	`\mu M\cdot \text{min}^{-1}`
`K_{m,GPI,f}`	309	`\mu M`
`V_{max,GPI,r}`	900000	`\mu M\cdot \text{min}^{-1}`
`K_{m,GPI,r}`	900	`\mu M`
`V_{max,PGM,f}`	7300	`\mu M\cdot \text{min}^{-1}`
`K_{m,PGM,f}`	170	`\mu M`
`V_{max,PGM,r}`	1000	`\mu M\cdot \text{min}^{-1}`
`K_{m,PGM,r}`	5000	`\mu M`
`K_{cat,UGPase,f}`	0.6	`s^{-1}`
`K_{S,UGPase,f}`	1000	`\mu M`
`K_{m1,UGPase,f}`	742	`\mu M`
`K_{m2,UGPase,f}`	42.4	`\mu M`
`K_{cat,UGPase,r}`	6.35	`s^{-1}`
`K_{S,UGPase,r}`	1000	`\mu M`
`K_{m1,UGPase,r}`	3040	`\mu M`
`K_{m2,UGPase,r}`	376	`\mu M`
`K_{cat,AGPase,f}`	0.6	`s^{-1}`
`K_{S,AGPase,f}`	1000	`\mu M`
`K_{m1,AGPase,f}`	742	`\mu M`
`K_{m2,AGPase,f}`	42.4	`\mu M`
`K_{cat,AGPase,r}`	6.35	`s^{-1}`
`K_{S,AGPase,r}`	1000	`\mu M`
`K_{m1,AGPase,r}`	3040	`\mu M`
`K_{m2,AGPase,r}`	376	`\mu M`
`K_{cat,SS}`	1.05	`s^{-1}`
`K_{m,SS}`	220	`\mu M`
`V_{max,SP}`	8530	`\mu M\cdot \text{min}^{-1}`
`K_{S,SP}`	1000	`\mu M`
`K_{m1,SP}`	360	`\mu M`
`K_{m2,SP}`	11940	`\mu M`
`V_{max,SPS}`	6000	`\mu M\cdot \text{min}^{-1}`
`K_{S,SPS}`	1000	`\mu M`
`K_{m1,SPS}`	900	`\mu M`
`K_{m2,SPS}`	2000	`\mu M`
`V_{max,SPP}`	15700	`\mu M\cdot \text{min}^{-1}`
`K_{m,SPP}`	350	`\mu M`

[1] Dilke Kazan, Agnes Camurdan, Amable Hortacsua. The effect of glucose concentration on the growth rate and some intracellular components of a recombinant E. coli culture. Process Biochemistry. 1995; 30(3):269-273. doi: 10.1016/0032-9592(95)85008-2.

[2] Broddrick, J. T., Rubin, B. E., Welkie, D. G., Du, N., Mih, N., Diamond, S., Lee, J. J., Golden, S. S., & Palsson, B. O. (2016). Unique attributes of cyanobacterial metabolism revealed by improved genome-scale metabolic modeling and essential gene analysis. Proceedings of the National Academy of Sciences of the United States of America, 113(51), E8344–E8353. https://doi.org/10.1073/pnas.1613446113

[3] Kim YJ, Park JY, Suh SH, Kim MG, Kwak HS, Kim SH, Heo EJ. Development and validation of a predictive model for pathogenic Escherichia coli in fresh-cut produce. Food Sci Nutr. 2021 Oct 29;9(12):6866-6872. doi: 10.1002/fsn3.2642. PMID: 34925814; PMCID: PMC8645737.

[4] Sabri S, Nielsen LK, Vickers CE. Molecular control of sucrose utilization in Escherichia coli W, an efficient sucrose-utilizing strain. Appl Environ Microbiol. 2013 Jan;79(2):478-87. doi: 10.1128/AEM.02544-12. Epub 2012 Nov 2. PMID: 23124236; PMCID: PMC3553775.

[5] Dilek Kazan, Agnes Camurdan, Amable Hortacsu. The Effect of Glucose Concentration on the Growth Rate and some Intracellular Components of a Recombinant E. coli Culture. Process Biochemistry Vol. 30, No. 3,pp. 269-273.1995

[6] Ropers D, Baldazzi V, de Jong H. Model reduction using piecewise-linear approximations preserves dynamic properties of the carbon starvation response in Escherichia coli. IEEE/ACM Trans Comput Biol Bioinform. 2011 Jan-Mar;8(1):166-81. doi: 10.1109/TCBB.2009.49. PMID: 21071805.

[7] Christoph Marschall, Regine Hengge-Aronis. Regulatory characteristics and promoter analysis of csiE, a stationary phase-inducible gene under the control of `\sigma`S and the cAMP-CRP complex in Escherichia coli. https://doi.org/10.1111/j.1365-2958.1995.mmi\_18010175.x

[8] Silva, I.J., Saramago, M., Dressaire, C., Domingues, S., Viegas, S.C. and Arraiano, C.M. (2011), Importance and key events of prokaryotic RNA decay: the ultimate fate of an RNA molecule. WIREs RNA, 2: 818-836. https://doi.org/10.1002/wrna.94

[9] Raychaudhuri, A., Jerga, A., & Tipton, P. A. (2005). Chemical mechanism and substrate specificity of RhlI, an acylhomoserine lactone synthase from Pseudomonas aeruginosa. Biochemistry, 44(8), 2974–2981. https://doi.org/10.1021/bi048005m

[10] Parsek, M. R., Val, D. L., Hanzelka, B. L., Cronan, J. E., Jr, & Greenberg, E. P. (1999). Acyl homoserine-lactone quorum-sensing signal generation. Proceedings of the National Academy of Sciences of the United States of America, 96(8), 4360–4365. https://doi.org/10.1073/pnas.96.8.4360

[11] Maeda, H., Fujita, N., & Ishihama, A. (2000). Competition among seven Escherichia coli sigma subunits: relative binding affinities to the core RNA polymerase. Nucleic acids research, 28(18), 3497–3503. https://doi.org/10.1093/nar/28.18.3497

[12] 张晓翠.(2018).人工构建的细菌群体感应系统动力学研究(硕士学位论文,厦门大学).

[13] van Zundert, G.C.P., Rodrigues, J.P.G.L.M., Trellet, M., Schmitz, C., Kastritis, P.L., Karaca, E., Melquiond, A.S.J., van Dijk, M., de Vries, S.J., and Bonvin, A.M.J.J. (2016). The HADDOCK2.2 Web Server: User-Friendly Integrative Modeling of Biomolecular Complexes. Journal of Molecular Biology 428, 720-725. https://doi.org/10.1016/j.jmb.2015.09.014.

[14] Kawano, F., Suzuki, H., Furuya, A., and Sato, M. (2015). Engineered pairs of distinct photoswitches for optogenetic control of cellular proteins. Nature Communications 6, 6256. 10.1038/ncomms7256.