Part I. Heat Transfer Model of Meal Replacement Machine
Motivation
The temperature has a complex effect on DNA replication, RNA transcription, synthesis of proteins, and other biochemical processes in the cell by affecting the activity of enzymes, the diffusion of substances, etc. Low temperatures reduce the activity of the engineered bacteria, while high temperatures may lead to the death of the engineered bacteria. Genetic circuits that respond to temperature have also been found in bacteria, which allow them to adapt to different temperature environments. In this project, we have designed engineered bacteria capable of responding independently to two temperature states and two light signals. The temperature-switching process of the culture fluid is very important for us to control the transition of the production state of the engineered bacteria. The lag time of the engineered bacteria in response to temperature and the residence time in different states together determine the product taste. The behavior of engineering bacteria during state switching is complex. The state of production and productivity is only easily measured and controlled if maintained in one state. Therefore, the uniform temperature distribution in the culture and the rapid temperature changeover are important for different flavors.
Introduction
We used the multi-physics field simulation software—Comsol—to model the meal replacement machine and solved the temperature variation during the heating process of the meal replacement machine by the finite element method. We investigated the effects of heating pad temperature, stir, and heat transfer on the heating rate of the culture solution.
Figure 1. Model of the meal replacement machine.
Model
The meal replacement machine (Figure 1) is considered to be a well-insulated cubic box (with a side length of 20 cm). The reactor is a conical flask with a height of 14 cm and a thickness of 1.3 mm. The reactor was filled with 3 cm high water to simulate the culture solution. The rest of the model is filled with air. Four thermostatic heating plates are mounted symmetrically on the outer wall of the conical flask. A space of 1 mm between the thermostatic heating sheet and the conical flask is used to simulate the thermal conductivity in different thermal contact situations. Thermal conductivity, density, etc. are set using the default settings in the software's material library except for the data mentioned.
To avoid discontinuous temperature changes leading to non-convergence of the results, the initial temperature of the heating plates is set to the ambient temperature and then rapidly ramps up to the target temperature in 3s. We use a magnetic stirrer to ensure a uniform temperature of the culture solution in the meal replacement machine. Simulation of the fluid flow in a stirring process involves fluid mechanics calculations of very high complexity. Considering the computational cost of the model, the homogenization of the temperature distribution due to stirring in the conical flask is equivalently replaced in our model by increasing the thermal conductivity of the fluid. The remaining parameters can be easily modified in the material parameters.
Results
We designed the following seven simulations to investigate the effects of "thermostatic heating plate temperature", "stirring" and " heat transfer coefficient of the heating plate-conical flask” on the heating rate of the culture solution.
Table 1. Parameter settings for simulations
| Number | Heater Temperature | Stirred | Heat Transfer Coefficient |
|---|---|---|---|
| 1 | 353 K | Y | 0.025 |
| 2 | 353 K | Y | 0.1 |
| 3 | 353 K | Y | 0.25 |
| 4 | 353 K | Y | 0.05 |
| 5 | 353 K | N | 0.05 |
| 6 | 343 K | Y | 0.05 |
| 7 | 333 K | Y | 0.05 |
Figure 2. Temperature rise curve of culture fluid when heated by heating plates of
different temperatures.
We simulated the heating of the culture solution with thermostatic heating plates at different temperatures (60°C, 70°C, 80°C). The simulation results are in line with our expectations: Higher temperature heating plates will result in faster heating rates. To achieve a better heating effect, we should use higher-temperature thermostatic heating plates. In practice, a higher temperature of the heating plate will raise the temperature of the reactor wall. The engineered bacteria will be killed by the high temperature when they touch the reactor wall. Therefore, the heating temperature should be increased taking into account the safety of the engineered bacteria.
Figure 3. Temperature rise curve of the culture solution with and without stirring.
To maintain a uniform temperature in all parts of the culture solution, we added a stirring module to the meal replacement machine. We used the intermittent stirring strategy to avoid the death of engineered bacteria caused by stirring. Stirring increases convection in the culture (allowing sufficient heat exchange between the various parts of the culture, avoiding overheating of the culture near the reactor walls) so that the culture can maintain a relatively uniform temperature at all times. Considering the heat transfer process from the reactor to the culture fluid at the boundary between the reactor and the culture fluid, a higher temperature gap can lead to more heat transfer under the same conditions (according to the heat transfer equation). Without stirring, the heat transfer to the inside of the culture is slow. Therefore, heat accumulates at the boundary (Figure 4), leading to a higher temperature at the boundary than under full stirring. The temperature gap between the heating plate and the boundary is small, which is not helpful for heat transfer to the culture solution. It takes more time to raise the average temperature of the culture solution to the set temperature without stirring. With full agitation at all times, it takes only about 700 s to raise the culture temperature from about 26 °C to the target temperature of about 37 °C. This process takes more than 1000 s without stirring. In the experimental measurements, this process took about 840 s. Because our meal replacement machine uses an intermittent stirring strategy, the actual time used for heating is between "full stirring" and "no stirring".
Figure 4. Temperature distribution on the longitudinal section with/without
stirring.
Figure 5. Temperature rise curve of the culture solution under different thermal
contact (equivalent thermal conductivity: J/(m K)) between the heating plates and the reactor.
In our meal replacement machine, we use a thermostatic heating plate to heat the culture solution. The thermal contact between the thermostatic heating plate and the glass reactor has a significant effect on the temperature rise rate of the culture solution. We simulated the heating process for four thermal conductivities (0.025 J/(m K), 0.05 J/(m K), 0.1 J/(m K), 0.25 J/(m K)). (Dashed lines in figure 5: darker colors correspond to larger thermal conductivities) The four thermal conductivities correspond to "no contact between the heater and the reactor (thermal conductivity close to air)", "partial thermal contact between the heater and the reactor (0.05 J/(m K), 0.1 J/(m K))" and "Good contact between heating sheet and reactor (thermal conductivity close to reactor wall)". The simulation results show that the thermal contact between the heating plate and the reactor of our current meal replacement machine is not good enough. By improving the thermal contact effect, we can theoretically increase the temperature rise speed by up to 4 times. Thus, the time required for temperature rise can be controlled within 3 min 30 s, which can significantly improve the user comfort and the accuracy of meal replacement taste control!
Part II. Model for the Bistable Switch
Variable definitions
Firstly, we provide the meaning of abbreviations.
Table 1. Definitions of parameters
| Parameter | Description |
|---|---|
| $M_l$ | The content of LaclLOV mRNA |
| $M_c$ | The content of cl:LVA mRNA |
| $M_y$ | The content of YFP mRNA |
| $M_m$ | The content of mcherry mRNA |
| $P_{l}$ | The content of LaclLOV |
| $P_{c}$ | The content of cl:LVA |
| $P_{y}$ | The content of YFP |
| $P_{m}$ | The content of mcherry |
| $k_{l}$ | Translation rate of LaclLOV |
| $k_{c}$ | Translation rate of cl:LVA |
| $k_{y}$ | Translation rate of YFP |
| $k_{m}$ | Translation rate of mcherry |
| $\alpha_{l}$ | Transcription rate of mRNA of LaclLOV |
| $\alpha_{c}$ | Transcription rate of cl:LVA |
| $\alpha_{y}$ | Transcription rate of YFP |
| $\alpha_{m}$ | Transcription rate of mcherry |
| $d_{l}$ | Degradation rate of LaclLOV |
| $d_{c}$ | Degradation rate of cl:LVA |
| $d_{y}$ | Degradation rate of YFP |
| $d_{m}$ | Degradation rate of mcherry |
| $\gamma_{l}$ | Degradation rate of LaclLOV mRNA |
| $\gamma_{c}$ | Degradation rate of cl:LVA mRNA |
| $\gamma_{y}$ | Degradation rate of YFP mRNA |
| $\gamma_{m}$ | Degradation rate of mcherry mRNA |
| $\beta_{ltocm}$ | Inhibition coefficients of LaclLOV on promoters of cl:LVA and mcherry |
| $\beta_{ctoy}$ | Inhibition coefficients of cl:LVA on promoter of YFP |
| $C_{a}$ | Content of transcribable DNA |
| $C_{d}$ | Content of blocked DNA |
| $C_{f}$ | Content of free state DNA |
| $k_{fa}$ | Positive reaction rate constant for the transition from the free state DNA to the transcribable state |
| $k_{fa}^{-1}$ | Reverse reaction rate constant for the transition from the free state DNA to the transcribable state |
| $k_{fd}$ | Positive reaction rate constant for the transition from the free state DNA to the non-transcribable state |
| $k_{fd}^{-1}$ | Reverse reaction rate constant for the transition from the free state DNA to the non-transcribable state |
| $P$ | Blocking protein content |
| $v_{tr}$ | True transcription rate |
| $v_{maxtr}$ | Maximum transcription rate |
Background of Lacl:LOV
For $mRNA$, we have
\begin{equation}\label{eq1}\frac{dM_l}{dt}=\alpha_{l}-\gamma_{l}M_l \end{equation}
It reach equilibrium quickly, so
\begin{equation}\label{eq2}\frac{dM_l}{dt}=0\Rightarrow M_l=\frac{\alpha_{l}}{\gamma_{l}}\end{equation}
We used Lacl:LOV novel light fusion protein and cI inhibitory factor with LVA degradation tag to constitute a bistable switch. Two different proteins are expressed in the absence and in the presence of blue light. In our system, Lacl:LOV is controlled by a strong promoter and has no factors that repress its expression. Based on the data from the experimental group, we only consider protein expression after the bacterial content is stabilized, so that its DNA content is considered as 1 and the protein content rapidly reaches equilibrium. Based on the underlying ODE, we derived that
\begin{equation}\label{eq3}\frac{dP_{l}}{dt}=k_{l}M_l-d_{l}P_{l}\end{equation}
So we bring (2) into (3) and get
\begin{equation}\label{eq4}\frac{dP_l}{dt}=\frac{k_l\alpha_l}{\gamma _{l}}-d_{l}P_l\end{equation}
We also assume that the reaction reaches equilibrium extremely quickly, so
\begin{equation}\label{eq5}\frac{dP_l}{dt}=0\Rightarrow P_l=\frac{k_l\alpha_l}{\gamma_l d_l}\end{equation}
Changes of DNA
Three types of DNA can be considered: free DNA that is neither transcribed nor repressed by the promoter, non-transcribed DNA that is repressed by transcription, transcriptional DNA that is combined with a helicase, the reaction from free DNA to transcriptional DNA and transcription of transcriptional DNA, all of which are zero-level reactions (sufficient for helicases, etc.), and free DNA that is repressed by a repressor to a non-transcribed state.
$P$ combined with $C_{f}$ to make it non-transcribable, and according to the previous derivation and the set parameters, we get the following conversion formula
\begin{equation}\label{eq6}C_{a}\stackrel{k_{fa}}\rightleftharpoons C_{f} \stackrel{k_{fd},nP}\rightleftharpoons C_{d}\end{equation}
The rate of the inverse reaction is shown in the corresponding parameter table, combined with our previously derived chemical equilibrium equation
\begin{equation}v_{tr}=\frac{v_{maxtr}k_{fa}}{k_{fa}+1+k_{fd}P^{n}}\end{equation}
We replaced some parameters$\;\;\;\;$
$\beta_{fd}=\frac{1+k_{fa}}{k_{fa}},V_{maxtr}=\frac{v_{maxtr}k_{fa}}{k_{fa}}$，we got
\begin{equation}v_{tr}=\frac{V_{maxtr}}{\beta_{fd}+P^{n}}\end{equation}
And our differential equation for $mRNA$ takes the form
\begin{equation}\frac{dM_{i}}{dt}=\frac{\alpha_{i}}{\beta_{fd_i}+P^{\beta}}-\gamma_{i}M_{i},\end{equation}
where $i=m$, $c$, $y$, respectively.
Notes: In the following part, we assumed that $\beta_{fd_i}\approx 1$ because $k_{fa}>>1$.
cl:LVA, mCherry, YFP
In the absence of blue light, clLVA promoter Ptrc-2 is repressed by LocI:LOV protein with coefficient $\beta_{ltocm}$. Denote that the maximum activity coefficient and degradation rate of the promoter is $\alpha_c$ and $\gamma _c$, respectively, then mRNA ($M_c$) production rate is
\begin{equation}\frac{dM_c}{dt}=\frac{\alpha_c}{1+P_l^{\beta_{ltocm}}}-\gamma_cM_c\end{equation}
According to ODE, denote the comtent of cl:LVA as $P_c$, the transcription rate as $k_c$, the degradation as $d_c$, so we have
\begin{equation}\frac{dP_c}{dt}=k_cM_c-d_cP_c\end{equation}
them combined the last step we can get the content of cl protein when reaching balance.
At equilibrium, LacI:LOV protein represses the upstream promoter Ptrc-2 in the absence of blue light so that it conbined with the repressor and could not be transcribed. With blue light the repression is lifted, so we considered $\beta$ to denote the repression coefficient, and when $\beta=0$, the repression is lifted.
\begin{equation}\frac{dM_m}{dt}=\frac{\alpha_m}{1+P_l^{\beta_{ltocm}}}-\gamma_mM_m\end{equation}
\begin{equation}\frac{dP_m}{dt}=k_mM_m-d_mP_m\end{equation}
cILVA can regulate the promoter based on the phage $\lambda$ pR promoter, and the PR promoter has two DNA binding sites for $\lambda$ cI blockers, whose activity is inhibited after binding.
\begin{equation}\frac{dM_y}{dt}=\frac{\alpha_y}{1+P_c^{\beta_{ctoy}}}-\gamma_yM_y\end{equation}
\begin{equation}\frac{dP_y}{dt}=k_yM_y-d_yP_y\end{equation}
Parameter simplification
We parametrically simplif the first set of equations for mRNA by bringing in the value of $P_l$ at equilibrium to obtain two first-order linear differential equations, and we let $\phi _l$ denote $\frac{k_l \alpha_l}{d_l \gamma_l}+1$ and solve the two mRNA rate equations
\begin{equation}M_c=A_c(e^{B_ct}-1),\end{equation}
\begin{equation}M_m=A_m(e^{B_mt}-1).\end{equation}
Bringing it into the ODE equation, we obtained the differential equation for the target product $P_m$, and we then performed a parametric simplification such that $R_m$ denotes $\frac{P_m}{A_c}$ , simplifying to obtain another set of differential equations
\begin{equation}\frac{dR_m}{dt}=e^{B_ct}-d_mR_m-1.\end{equation}
So we can get the rate equation of $P_c$
\begin{equation}\frac{dR_c}{dt}=e^{B_ct}-d_cR_c-1.\end{equation}
Results
We modeled each nutrient of the three plasmids and the YFP-tagged material in the same way and obtained the production rates of the other three nutrients
\begin{equation}\frac{dR_x}{dt}=e^{B_xt}-d_xR_x-1\end{equation}
using MATLAB, we solved the differential equations and got the results
\begin{equation}R_{x}=\frac{B_{x}e^{-d_xt}}{d_x^2+AB}-\frac{e^{d_xt+B_xt}}{d_x+B_x}-\frac{1}{d_x}.\end{equation}
Simulation
We built a simbiology model to simulate the four states and switches. To demonstrate the effect of the presence or absence of blue light on the target product content, we used YFP fluorescence to show the non-productive (called OFF hereinafter, and the other state is called ON) state. In addition, we considered that the inhibition coefficient is 1 in the presence of blue light and 0 in the absence of blue light, and the Lacl:LOV content is constant, so the inhibition part $1+lacl:LOV^\beta_ltocm$ is a constant, so we simplified the parameters so that $k_{fc}$ replaced $\frac{\alpha_c}{1+Lacl:LOV^\beta_{ltocm}}$ to indicate the effect of blue light.
We set the temperature as a constant, and our bacteria were initially set to the OFF state. In our simulation, we successively switched the state to ON and OFF at 50min and 500min, and our results are as follows.
Figure 1. Overview of system behavior under the control of bistable switch.
Figure 2. System behavior after switching to ON state.
Figure 3. System behavior after switching to OFF state.
For mRNA concentration, we changed the states of switches quickly, and the result showed that mRNA equilibrium can be reached within 10 minutes for both productive and non-productive states.
Figure 4. Changes of [mRNA] when state switches.
However, the result above also implied that when bacteria are in the OFF state, mCherry might produce slowly, and we call this effect “leakage”. To evaluate the intensity of such “leakage”, numerical simulation was used to calculate the concentration of mCherry at the OFF state, which is shown below.
Figure 5. “Leakage” at non-productive state.
It showed that if we store our bacteria at room temperature, mCherry will reach a relatively low concentration. However, as the production process usually lasts for not longer than two days, such “leakage” can be neglected
.Optimization
To simplify the notations, we denote that for component $i$, its production speed can be treated as a constant $v_{pi}$ (as mRNA equilibrium can be reached in a relatively short time, see Figure 4) and degrade rate is $k_{di}$ (we assume that it is unrelated to temperature), $i=1,2,3$. Denote that the concentration of component $i$ is $c_i$, then
\begin{equation}\dfrac{dc_i}{dt}=\begin{cases}v_{pi}-k_{di}c_i, & \text{state=i}, \\ -k_{di}c_i, & \text{otherwise}.\end{cases}\end{equation}
According to our experiment, transition time between two states cannot be neglected, therefore, we use $\tau_{ij}$ to describe time consumption from state $i$ to state $j$, and assume that there is no production when the state switches. The goal of production is to let $c_i=c_{fi}$ eventually. We firstly consider the order $1\rightarrow 2\rightarrow 3$. We assume time consumption at each state is $t_i$, by solving these ODEs, we derive
\begin{align}t_3 &= \dfrac{\ln v_{p3} - \ln(v_{p3}-c_{f3}k_{d3})}{k_{d3}}, \\ t_2 &= \dfrac{\ln v_{p2} - \ln(v_{p2}-c_{f2}k_{d2}e^{k_{d2}(t_3+\tau_{23})})}{k_{d2}}, \\ t_1 &= \dfrac{\ln v_{p1} - \ln(v_{p1}-c_{f1}k_{d1}e^{k_{d1}(t_2 + t_3 + \tau_{12} +\tau_{23})})} {k_{d1}}.\end{align}
Therefore, we can estimate the cost of this production process, denoted as $u_{i,j,k}$, where $i,j,k$ is some permutation of $1,2,3$. For example, if we estimate total time consumption, we have
\begin{equation}u_{i,j,k}^t=t_1+t_2+t_3+\tau_{12}+\tau_{23}.\end{equation}
Alternatively, if we estimate total energy consumption, and we denote $P_i$ as power state $i$, and $W_{ij}$ as energy consumption when switching from state $i$ to state $j$, then
\begin{equation}u_{i,j,k}^W=P_1t_1+P_2t_2+P_3t_3+W_{12}+W_{23}.\end{equation}
By estimating $u_{i,j,k}$ for all 6 permutations, our project can choose the best order to produce the required meals.
Part III. Shortest Time model
Consider a system described by the equation below:
\begin{equation}\dot{x_i}=\begin{cases}k_i-c_ix_i, & T_{i-1}\leq t < T_i, \\ -c_ix_i, & \text{otherwise},\end{cases}\qquad T_i=\sum_{j=1}^it_i\qquad x_i(0)=0\qquad i=1,2,3,4.\end{equation}
Now we need to determine $t_i$ such that $x_i(T_4)=f_i$. A naive solution is as follows. As
\begin{equation}x_4=\dfrac{k_4}{c_4}(1-\mathrm{e}^{-c_4(t-T_3)}),\end{equation}
we can derive that
\begin{equation}t_4=\dfrac{\ln k_4-\ln(k_4-f_4c_4)}{c_4}.\end{equation}
For $t_3$, if $x_3(T_4)=f_3$, then $x_3(T_3)=f_3\mathrm{e}^{c_3t_4}$, therefore,
\begin{equation}t_3=\dfrac{\ln k_3-\ln(k_3-f_3c_3\mathrm{e}^{c_3t_4})}{c_3}.\end{equation}
And that is similar for $t_2$ and $t_1$, that is,
\begin{equation}t_2=\dfrac{\ln k_2-\ln(k_2-f_2c_2\mathrm{e}^{c_2(t_3+t_4)})}{c_2},\end{equation}
\begin{equation}t_1=\dfrac{\ln k_1-\ln(k_1-f_1c_1\mathrm{e}^{c_1(t_2+t_3+t_4)})}{c_1}.\end{equation}
Naively, we only need to compare the results of 24 permutations to derive the shortest time.
Acknowledgements
