Model

1.Abstract

The substances sensed by the biosensors EtnR1 and EtnR2 are ethylene oxide, and they cannot directly sense ethylene. Therefore, in order for the biosensor of this project to realize ethylene detection, it is necessary to introduce and express an enzyme that converts ethylene into ethylene oxide in a different strain of E. coli from that which expresses EtnR1 and EtnR2. However, since the enzyme used in this study originates from Burkholderia cepacia G4 and its expression level in E. coli is unknown, we constructed an enzyme expression prediction program using Python. In addition, the ethylene monooxygenase used in this study is a mutant of toluene o-xylene monooxygenase derived from the above bacterium. Because prior studies on this enzyme are scarce, we conducted comparisons of mutant subunits using UCSF Chimera, as well as predictions of epoxidation reactions during enzyme expression, in order to clarify the structural mutations of the enzyme.

2.Introduction

In developing the present prediction program, it was found that the V106F and A113F mutants of toluene/o-xylene monooxygenase (Tom) can essentially function as ethylene monooxygenases(1).Therefore, to investigate the effect of amino acid mutations on the three-dimensional structure of the enzyme, the mutant and wild-type subunits were compared using UCSF Chimera; however, due to time constraints, the evaluation of enzyme–substrate binding affinity and docking scores through molecular docking could not be conducted.

3.Prediction of Gene Expression and Structural Comparison of Mutant Subunits

3-1. Prediction of Enzyme Expression Levels

3-1-1. Preparation of Machine Learning Data

To construct the expression prediction program in Python, it was necessary to prepare training data consisting of nucleotide sequences with known expression levels in E. coli. Therefore, we searched Escherichia coli str. K-12 substr. MG1655 in PaxDb (2), and obtained protein names and expression levels from the “E. coli – Whole organism (Integrated)” abundance list. Next, we searched the corresponding proteins in UniProt (3), and from the search results with the organism Escherichia coli (strain K12), we followed the GenBank links to obtain DNA nucleotide sequences for the training dataset.

3-1-2. Construction of Python Code

For the DNA sequences and expression levels prepared in 3-1-1, we extracted MFE and secondary structures using RNAfold, and transformed them into 768-dimensional vectors with Longformer to construct features (MFE + vector). These were divided into 80% training and 20% test sets, and a random forest regression model was trained to predict expression levels. Furthermore, predictions were performed for unknown sequences (TomA0–TomA5) as well as known sequences (the nucleotide sequences of rplF and rpsP), and we created bar graphs of feature importance and scatter plots of observed vs. predicted values to evaluate model performance and the contribution of each feature.

3-1-3. in silico Prediction Results of Expression Levels

The prediction results of expression levels are shown below (Tables 1 and 2), together with the bar graph of feature importance (Graph 1) and the scatter plot of observed vs. predicted values (Graph 2). As shown in Tables 1 and 2, when the training data consisted of only 40 entries, the test MSE was large and the accuracy was low. Therefore, by newly adding 57 pairs of sequences and expression levels, and further removing the two pairs with the highest expression as outliers, we obtained the results shown in Tables 3 and 4. In addition, changes were observed in the bar graphs of feature importance and the scatter plots of observed vs. predicted values (Graphs 3 and 4).

From the graphs, it can be seen that when the training data consisted of only 40 samples, the MSE was extremely large, and there was a significant difference between the observed and predicted values even for known sequences. The main reason is likely that the dimensionality of each feature exceeded the number of samples, which caused the prediction model to memorize noise or accidental patterns in the data, resulting in overfitting. As an improvement, we increased the number of samples and removed the two pairs with the highest expression levels as outliers. As a result, the Test MSE decreased and the prediction accuracy improved. Furthermore, compared to the case with 40 samples, the correlation between observed and predicted values appeared to be stronger.

3-2. Structural Comparison among Subunits

3-2-1. Structural Comparison among Subunits

First, as reference information, we obtained structural data for toluene/o-xylene monooxygenase from different species (4) as well as structural data of ethylene (5). Ideally, it would have been preferable to use the structural data of Burkholderia cepacia G4, but since such information was not available, we referred to enzymes from other species. Nevertheless, this enzyme is also composed of three subunits, similar to the mutant enzyme in this study, and was therefore considered useful for understanding the overall structure of (αβγ) (Figure 1). Subsequently, we attempted to predict the tertiary structure of the target enzyme using AlphaFold2 (6) on Google Colab. However, when we tried to run the prediction for the (αβγ)₂ complex, the process was stopped due to insufficient GPU resources. Therefore, we limited the prediction to TomA3, which contains the catalytic active site.

3-2-2. Structural Comparison of TomA3 Subunits Using UCSF Chimera

Using AlphaFold2, we predicted the tertiary structures of TomA3 (monomer) of A113F, V106F, and Wildtype, and obtained PDB files (Figures 3–5). In these predictions, AlphaFold2 stored the pLDDT values in the B-factor column, and the reliability was visualized by color coding. Therefore, when displaying the output files in UCSF Chimera (7), we applied the B-factor command to show the reliability in the same manner (Table 5). Subsequently, we displayed the mutant TomA3 subunits together with the existing structure of ethylene in UCSF Chimera (Figures 6 and 7).

From these figures, based on the comparison of residues 106 and 113 with markers (Figures 3–5), we found that although no remarkable structural differences were observed visually, the predicted structural reliability based on pLDDT was higher for Wildtype and V106F near residue 113 than for A113F, and higher for Wildtype than for V106F near residue 106.

This may be because the mutations examined here are not well known, and thus there might have been discrepancies between AlphaFold2’s training data and the actual structures. Another possibility is that the side chain of phenylalanine (Phe/F) has a greater molecular mass than those of alanine (Ala/A) or valine (Val/V), potentially causing steric hindrance by colliding with neighboring amino acid residues. Alternatively, when considered together with previous literature (8), it is suggested that the main driving force for polypeptide folding in aqueous environments is hydrophobic interaction. Since phenylalanine has a large aromatic ring side chain, it can form stronger hydrophobic interactions than alanine or valine, leading to a smaller active site. As a result, it is possible that the enzyme became able to catalyze monooxygenase reactions not only for large molecules such as toluene and o-xylene, but also for smaller molecules such as ethylene.

Previous studies have reported that the reaction rate of ethylene monooxygenase differs greatly between V106F and A113F (1). However, due to time constraints, we were unable to evaluate binding affinity or docking scores through molecular docking in this study. Therefore, in future work, we aim to examine these factors, identify the determinants of catalytic rates, and link them to wet-lab experiments.

Figure 6: Active site and ethylene in V106F

Figure 7: Active site and ethylene in A113F

4. Prediction of Epoxidation Reaction during Enzyme Expression

4-1. Reaction Mechanism

The reaction mechanism of TOM is shown in Figure 8 (1). As can be seen in Figure 8, TOM possesses a reaction mechanism that epoxidizes aromatic compounds. However, in the case of its variants, V106F and A113F, the only difference lies in the fact that the substrate is replaced with ethylene, and it is expected that the same reaction mechanism as the normal TOM occurs (1). Nevertheless, in the reaction mechanism of TOM, chemical equations describing the reaction, as shown in Figure 8, have not been fully clarified in existing studies. On the other hand, for alkene monooxygenases that oxidize alkenes, including ethylene, the detailed chemical equations have already been elucidated. Therefore, it is assumed that a similar reaction mechanism occurs in V106F and A113F. Consequently, the reaction mechanism in Figure 8 can be expressed by the following chemical equation1 (11).

equation1 $$ \text{NAD(P)H} + C_2H_4 + O_2 + H^+ \rightarrow \text{NAD(P)}^+ + C_2H_4O + H_2O $$

4-2. OED and Parameter

4-2-1. Parameter

The parameters appearing later are first defined here. For simplicity,

$$ A = [\mathrm{NADH}] \hspace{20pt} B = [\mathrm{O_2}] \hspace{20pt} C = [\mathrm{C_2H_4}] $$

the parameters are summarized in Table 7.

4-2-2. ODE

Based on the fact that the reaction in equation2. Is an ordered reaction (12), the three substrates \( [\mathrm{NAD(P)H}] \), \( [\mathrm{C_2H_4}] \), and \( [\mathrm{O_2}] \) were taken as parameters, while the concentration of the product \( [\mathrm{C_2H_4O}] \) was expressed as \( [\mathrm{C_2H_4O}] \). When the rate of change in \( [\mathrm{C_2H_4O}] \) is expressed as \( v \), the OED for \( v \) is described as follows:

equation2 $$ v = \dfrac{V_{\max} [A][B][C]} {K_{ia}K_bK_c + K_bK_c[A] + K_{ia}K_c[B] + K_{ia}K_b[C] + K_c[A][B] + K_b[A][C] + K_{bc}[B][C] + [A][B][C]} \ $$

Here, assuming that \( \mathrm{NAD(P)H} \) and \( \mathrm{O_2} \) exist in excess amounts, \( [\mathrm{NAD(P)H}] \) and \( [\mathrm{O_2}] \) are considered constant before and after the reaction. Therefore, equation2 can be approximated as the following first-order reaction equation for \( v \):

equation3 $$ \frac{d[\mathrm{C_2H_4O}]}{dt} = v([C]) = \frac{V_{\max}^\ast [C]}{K_m^\ast + [C]} \ $$

In addition, the rate of decrease of \( [\mathrm{C_2H_4}] \) is expressed as \( -v \), and similarly represented as follows:

equation4 $$ \frac{d[\mathrm{C_2H_4}]}{dt} = -v([C]) = -\frac{V_{\max}^\ast [C]}{K_m^\ast + [C]} \ $$

4-3. In silico Simulation (Initial Data)

4-3-1. Initial Values

The initial values of \( [\mathrm{V_{max}}] \), \( [\mathrm{K_m}] \), and \( [\mathrm{C}] \) in equation3 were set as shown below (Tables 8 and 9) (1).

4-3-2. Simulation

Based on the initial values set in [1] and the OED above, the reaction time was first set to 250 minutes, and the simulation was carried out (Graphs5 and 6).

[1] Graph5: Concentration changes of C₂H₄ and C₂H₄O catalyzed by the expressed A113F enzyme [2] Graph6: Concentration changes of C₂H₄ and C₂H₄O catalyzed by the expressed V106F enzyme

At this point, the graph in [2] Graph6 appeared linear and did not show the curve expected from a first-order reaction equation as in [1]. The reason for this result is presumed to be that the maximum reaction velocity of V106F is far smaller compared to A113F. Therefore, only for V106F, the reaction time was extended up to 10,000 minutes, and the graph was re-plotted. The results are shown below (Graph7).

Graph7: Concentration changes of C₂H₄ and C₂H₄O during the expression of V106F (reaction time extended to 10,000 minutes) <>It can be observed from the plotting results that a curve similar to that of A113F was obtained.

1. Abstract

In this project, we developed and analyzed a mathematical model describing the behavior of a co-culture system composed of Escherichia coli engineered to express monooxygenase, which converts ethylene into ethylene oxide, and sensor E. coli responsive to ethylene oxide. First, using single-strain cultivation data, we estimated fundamental parameters such as the maximum specific growth rate and the substrate half-saturation constant. Next, enzyme expression levels and catalytic turnover numbers were estimated, and the epoxidation reaction was described using the Michaelis–Menten framework. Furthermore, growth inhibition by the product, ethylene oxide, was incorporated into the model. Substrate consumption, product accumulation, and sensor response under co-culture conditions were then represented as a system of ordinary differential equations (ODEs). Through this model, we simulated the dynamics under different initial biomass and substrate concentrations and identified conditions that maximize the sensor response.

2. Introduction

In this project, the biosensor is designed on the premise of co-culturing two genetically transformed Escherichia coli strains with distinct functional roles. One strain expresses toluene o-xylene monooxygenase (Tom), which oxidizes ethylene into ethylene oxide (hereafter referred to as the Tom-expressing strain). The other strain expresses EtnR1 and EtnR2, thereby promoting transcription from the EtnP promoter in the presence of ethylene oxide (hereafter referred to as the sensor strain).

We modeled the growth dynamics of both strains using the Monod equation, while also constructing a system of ordinary differential equations (ODEs) to describe the induction of Tom expression by ethylene, the changes in ethylene oxide concentration, and the response of EtnP as a function of ethylene oxide concentration.

Furthermore, by applying this model, we simulated the magnitude of the EtnP response under varying initial cell density ratios. We then searched for the optimal set of initial conditions—[\(X_{0,1}, X_{0,2}, S_{0}\)], where \(X_{0,1}\) denotes the initial biomass of the Tom-expressing strain, \(X_{0,2}\) the initial biomass of the sensor strain, and \(S_{0}\) the initial substrate concentration—under which the EtnP response reaches its maximum.

3. Estimation of Fundamental Parameters through Single-strain Cultivation Using the Monod Equation

Simulation Results of the Tom-expressing Strain

Symbol	Value	Unit
\(\mu_{max}\)	0.7007	h⁻¹
\(K_S\)	0.7850	g·L⁻¹
\(Y_{X/S}\)	1.361×10⁹	cells·mL⁻¹·(g·L⁻¹)⁻¹

Simulation Results of the Sensor Strain

Symbol	Value	Unit
\(\mu_{max}\)	0.5439	h⁻¹
\(K_S\)	0.3220	g·L⁻¹
\(Y_{X/S}\)	1.702×10⁹	cells·mL⁻¹·(g·L⁻¹)⁻¹

4. Co-culture Simulation (Before Ethylene Addition)

Objective

In this chapter, we aim to construct a growth model for a single strain based on the Monod equation and to estimate fundamental parameters—namely, the maximum specific growth rate (\(\mu_{max}\)), the substrate half-saturation constant (\(K_s\)), and the yield coefficient (\(Y_{X/S}\))—from experimental data. These parameters provide essential values that form the foundation for subsequent co-culture simulations.

$$ \mu_i(S) = \mu_{i,\max},\frac{S}{K_{S,i}+S} $$

$$ \frac{dX_i}{dt} = \mu_i(S)X_i,\quad i=1,2 $$

$$ \frac{dS}{dt} = -\sum_{i=1}^{2}\frac{1}{Y_{X/S,i}}\mu_i(S)X_i $$

The system of ODEs derived from the Monod model was numerically integrated using SciPy’s solve_ivp function. Vectorized operations were implemented with NumPy, and growth curves as well as substrate consumption curves were visualized using Matplotlib. The initial conditions were set to \(S_{0}=1.6\) g/L and \(X_{0,1}=X_{0,2}=4.0 \times 10^{8}\) cells/mL. Since this analysis assumed conditions before ethylene addition, product formation and inhibition terms were not included.

Simulation Results

Discussion

The bacterial concentrations began to diverge gradually around 200 minutes and reached a stable difference at approximately 300 minutes.

Under co-culture conditions, the substrate is consumed by both strains, making the choice of initial conditions a major determinant of the final cell densities. Due to the difference in growth rates between the Tom-expressing strain and the sensor strain, extended cultivation resulted in a noticeable disparity in cell concentrations between the two populations.

Moreover, adjusting the initial ratio of cell densities (\(X_{0,1}:X_{0,2}\)) may enable control of the balance between the two strains at the time of ethylene addition. For example, increasing the initial density of the Tom-expressing strain could offset the growth advantage of the sensor strain, thereby creating a more balanced system.

In addition, the initial substrate concentration (\(S_{0}\)) also plays a critical role. If \(S_{0}\) is too low, sufficient biomass cannot be achieved; conversely, if \(S_{0}\) is too high, competition between the strains may become excessively intense. The next chapter will examine in detail how these initial conditions influence the overall system response after ethylene addition.

5. Co-culture after Ethylene Addition

Objective

In this chapter, we aim to model and simulate the behavior of the co-culture system following ethylene addition. Specifically, an ODE model was constructed that integrates three processes: (i) the conversion of ethylene into ethylene oxide by the Tom-expressing strain, (ii) the growth-inhibitory effects of ethylene oxide, and (iii) the sensing and response of the sensor strain to ethylene oxide. This integrated framework allows us to predict the overall performance of the biosensor system and provides a foundation for exploring optimal operating conditions.

5-1. Estimation of Monooxygenase Expression and \(k_{cat}\)

$$[E] = \frac{C_{crop}Xf}{MW}$$

$$k_{\mathrm{cat}} = \frac{V_{\max}}{[E]}(\frac{10^{-3}M}{mM}\frac{1s}{60min})$$

Symbol	Definition	Unit
\(X\)	Cell concentration	cells·mL⁻¹
\(C_{crop}\)	Soluble protein content per cell	mg·mL⁻¹·(cell·mL⁻¹)⁻¹
\(MW\)	Enzyme molecular weight	g·mol⁻¹
\(f\)	Mass fraction of enzyme within total protein	-
\([E]\)	Enzyme concentration	mol·L⁻¹
\(V_{max}\)	Maximum rate (experimental)	mM·min⁻¹
\(k_{cat}\)	Turnover number	s⁻¹

Parameters and Rationale

Parameter	Value	Basis
\(C_{crop}\)	0.22/8.0×10⁸ mg·mL⁻¹·(cell·mL⁻¹)⁻¹	Literature: 0.22 mg/mL/OD; 1 OD = 8×10⁸ cells/mL(3)
\(MW\)	223100 g·mol⁻¹	Calculated from sequence; subunit composition α₂β₂γ₂(2)
\(f\)	1.0	Upper estimate (assumes all protein is active enzyme)
\(V_{max,A}\)	0.019 mM/min	Table 2 (A113F)(4)
\(V_{max,V}\)	0.0007 mM/min	Table 3 (V106F)(4)

The total protein content per cell was calculated using the known value of 0.22 mg/mL/OD (corresponding to 8×10⁸ cells/mL/OD), and assuming \(f=1\), the enzyme concentration \([E]\) was estimated. By applying the molecular weight (MW = 231.1 kDa), \([E]\) was converted into mol/L units. The two distinct \(V_{max}\) values correspond to different Tom mutants, and \(k_{cat}\) was estimated for each case.

Results: Estimated Enzyme Concentration and \(k_{cat}\)

cells·mL⁻¹	\([E]\)	\(k_{cat,A}\)	\(k_{cat,V}\)
1.0×10⁹	1.23×10⁻⁶	0.2569	0.0095

Discussion

The estimated enzyme concentrations and \(k_{cat}\) values were indirectly derived from literature data. To obtain more accurate values, it is necessary to experimentally determine what proportion of expressed Tom proteins are catalytically active. Since the current estimation assumes \(f=1\) (all expressed Tom proteins are active), refining this assumption would yield more precise results.

5-2. Epoxyethane Production Rate Dependent on Monooxygenase Concentration (Michaelis–Menten Model)

$$\frac{d[\mathrm{C_2H_4O}]}{dt}= v(C)= \frac{V_{\max}^\ast \, C}{K_m^\ast + C}$$

$$\frac{d[\mathrm{C_2H_4}]}{dt}= -v(C)= -\frac{V_{\max}^\ast \, C}{K_m^\ast + C}$$

$$V_{max}=k_{cat}[E]$$

For A113F mutant

Parameter	Initial value
\(k_{cat}\)	0.2569
\(K_m\)	0.5 (mM)
\(C (= [C₂H₄])\)	0.7 (mM)
\(P (= [C₂H₄O])\)	0.0 (mM)

For V106F mutant

Parameter	Initial value
\(k_{cat}\)	0.0095
\(K_m\)	0.5 (mM)
\(C (= [C₂H₄])\)	0.7 (mM)
\(P (= [C₂H₄O])\)	0.0 (mM)

This simulation was based on predictions of epoxygenation reactions during Tom expression. The reaction rate from C to P was expressed by the Michaelis–Menten equation, with \(V_{max}\) varying according to enzyme concentration. Incorporating these considerations, \([E]\) was calculated and an enzyme concentration-dependent Michaelis–Menten model was formulated. All Tom variants used henceforth refer to the A113F mutant. The Monod parameters determined in Chapter 1 were adopted, with cultivation conditions set to \(S_0=1.6\) g/L and OD = 0.5.

Discussion

In this simulation, the production rate of epoxyethane was calculated from the time-dependent enzyme concentration of the Tom-expressing strain using the Michaelis–Menten framework. The results demonstrated that as cell density increased, \([E]\) rose proportionally, leading to an increase in \(V_{max}\) and thereby causing a marked acceleration of epoxyethane production during later cultivation stages. Consequently, under constant ethylene supply, epoxyethane production is expected to accelerate over time until limited by substrate supply.

However, the current model contains several simplifications. First, enzyme expression was assumed to scale instantaneously and proportionally with cell density, without accounting for transcriptional and translational delays or induction times. Second, transport, volatilization, and diffusion losses of ethylene and products were neglected, although these may constrain production rates in real systems. Furthermore, neither the cytotoxicity of epoxyethane nor its potential inhibitory effects on enzymes were considered, meaning reduced production rates at high concentrations could not be reproduced.

Future improvements should incorporate dynamic expression models including induction delays, gas–liquid transfer processes, and product inhibition to better capture realistic behavior. Additionally, experimental time-course data of ethylene and epoxyethane concentrations would refine parameterization and enhance predictive accuracy.

5-3. Growth Inhibition by Epoxyethane

In this section, the parameter \(K_P\), representing the concentration-dependent growth inhibition by epoxyethane (P), was estimated.

Two biomass growth datasets were prepared under different initial P concentrations, and the following Monod-based model with an inhibition term was applied:

$$\mu(S,[C_2H_4O])=\mu_{max}\frac{S}{K_S+S}\frac{1}{1+\frac{P}{K_{P}}}​$$

This system of ODEs describing the time evolution of \(X\) (biomass) and \(S\) (substrate) was integrated numerically using scipy.integrate.solve_ivp. Estimation of \(K_P\) was performed with scipy.optimize.curve_fit by minimizing the deviation between observed \(X(t)\) values and model predictions. Both time \(t\) and P concentrations were simultaneously used as explanatory variables, allowing joint optimization across the two conditions to achieve a more stable estimate.

The experiments employed the sensor strain, with the same initial substrate concentration across all conditions.

Experimental data:

Download Kp_fitting_EtnR1R2_OD_EtoON.csv

Parameter Estimation Result

$$K_P = 2.251$$

Discussion

This analysis successfully estimated the parameter \(K_P\), which characterizes growth inhibition as a function of epoxyethane concentration.

Nevertheless, several limitations remain. The dataset was restricted to a narrow range of initial concentrations, and broader coverage would improve the reliability of the estimation. Moreover, inhibition was modeled by a simple non-competitive function \(\frac{1}{1+P/K_{P}}\), while in reality, inhibition patterns may be more complex, potentially requiring models incorporating Hill coefficients or toxicity thresholds. Additionally, P concentration was assumed constant; however, in practice it fluctuates due to production, degradation, diffusion, and volatilization. Incorporating dynamic P levels would allow the model to better represent actual system behavior.

Future studies should re-estimate parameters using multiple concentration conditions and high-resolution growth data, while also testing different inhibition models. Such refinements are expected to enhance the predictive accuracy of the model.

6. Prediction of EtnP Expression

To predict the expression level from the EtnP promoter in the co-culture system after ethylene addition, we integrated the previously constructed models: the growth model, the enzyme expression model, the product formation model, and the inhibition model.

A Hill equation describing the response magnitude of EtnP as a function of ethylene oxide concentration, reported by Claudia F. Moratt et al., was incorporated as the final output of the model(5).

$$y=\frac{2857P^{1.22}}{P^{1.22}+1.96^{1.22}}$$

The sensor output, denoted as \(y'\), was defined as the product of this Hill function and the cell density of the ethylene oxide sensor strain, normalized to OD units:

$$y'=\frac{X_2}{8\times10^{8}}y$$

All equations were unified into a single ODE system and simulated accordingly.

Before Ethylene Addition (0–600 min)

$$\mu_i(S) = \mu_{i,\max}\,\frac{S}{K_{S,i}+S}​$$

$$\frac{dX_i}{dt} = \mu_i(S)X_i,\quad i=1,2$$

$$\frac{dS}{dt} = -\sum_{i=1}^{2}\frac{1}{Y_{X/S,i}}\mu_i(S)X_i$$

$$[E] = \frac{C_{crop}X_1f}{MW}$$

After Ethylene Addition (600–1000 min)

$$\mu_i(S,P) = \mu_{i,\max}\,\frac{S}{K_{S,i}+S}\frac{1}{1+\frac{P}{K_{P}}}\quad i=1,2$$

$$\frac{dX_i}{dt} = \mu_i(S)X_i,\quad i=1,2$$

$$\frac{dS}{dt} = -\sum_{i=1}^{2}\frac{1}{Y_{X/S,i}}\mu_i(S)X_i$$

$$ [E] = \frac{C_{crop} X_1f }{ MW }$$

$$\frac{dP}{dt}= \frac{k_{cat}[E] \, C}{K_m^\ast + C}$$

$$y'=\frac{X_2}{8\times10^{8}}\frac{2857P^{1.22}}{P^{1.22}+1.96^{1.22}}$$

The initial conditions were set to \(S_0 = 5.0\) g/L and \(X_{0,1} = X_{0,2} = 4.0\times10^{8}\) cells/mL. At 600 min, ethylene was added at a concentration of 0.7 mM.

Result

7. Exploration of Optimal Initial Conditions

In this analysis, the sensor response \(y'(t)\) in the co-culture after ethylene addition was evaluated using two criteria:

• The maximum value of \(y’\) between 600–1000 min
• The earliest time at which \(y’=10000\) was reached

The ODE system included the estimated Monod parameters, the enzyme concentration model, the Michaelis–Menten reaction rate, and the inhibition term. Using SciPy’s solve_ivp, the system was integrated over two consecutive phases: 0–600 min (without ethylene) and 600–1000 min (with ethylene).

Search Strategy

The search was conducted in two stages:

1. Coarse Grid Search
   • The initial substrate concentration \(S_0\) was varied from 1.0–10.0 g/L.
   • The constraint \(X_{0,1}+X_{0,2}=8.0\times10^9\) cells/mL was imposed, with \(X_{0,1}, X_{0,2}\) distributed in the range [\(0, 8.0\times10^9\)].
   • For each condition, both the maximum \(y'\) during 600–1000 min and the time to reach \(y'=10000\) were computed.
   • Top candidates were extracted as seeds.
2. Bayesian Optimization (Optuna + TPE Sampler)
   • The best grid-search conditions were used as seeds.
   • Optuna was employed to efficiently sample the parameter space, improving accuracy while avoiding local optima.

Discussion

The prediction of EtnP expression showed that, under standard initial conditions, the model successfully reproduced the delayed rise of sensor response \(y’\), following the increase in ethylene oxide concentration \(P\) due to elevated enzyme concentration \([E]\). This confirmed that the model consistently captures the timing and amplification of the response to ethylene addition.

By exploring initial conditions, combinations of \(S_0, X_{0,1}, X_{0,2}\) were identified that produced stronger and earlier responses compared with the standard conditions. This suggests that tuning the initial substrate concentration and cell ratio allows control over both sensitivity and response speed of the sensor system. Increasing the proportion of the Tom-expressing strain shortened the response time after ethylene addition, while increasing the proportion of the sensor strain amplified the response magnitude.

Nevertheless, several simplifications were present in the current model. Enzyme expression was assumed to be instantaneously proportional to cell density, neglecting transcriptional and translational delays as well as induction dynamics. Gas–liquid partitioning, diffusion losses, and metabolic degradation of ethylene and epoxyethane were not considered, which may cause actual responses at high concentrations to be smaller than predicted. Moreover, sensor response was modeled solely by a Hill equation, whereas stochastic models that incorporate transcription factor dynamics and biological noise may be required.

Future improvements should incorporate induction delays, mass transfer processes, and cytotoxic effects of products. Additionally, parameter refinement using experimental data will be necessary to enhance quantitative predictive accuracy.

1. Abstract

The concept of our biosensor centers on the prevention of foodborne illness caused by bacteria attached to fresh produce, the detection of over-ripening, and the visualization of these changes. In this section, “Evaluation of the Antimicrobial Activity of Nisin,” we assess the functionality of nisin, which plays the role of protecting produce from bacterial contamination within the biosensor system.

We constructed a multivariate model describing the antimicrobial activity of nisin as a function of three parameters: concentration, temperature, and pH. First, we estimated the relative impact of each parameter on the bactericidal efficacy of nisin. Subsequently, we evaluated its antimicrobial performance under the conditions anticipated in this project.

2. Introduction

We aimed to quantitatively elucidate how the bactericidal efficacy of the natural antimicrobial peptide nisin varies as a function of three key factors—concentration, pH, and temperature—and to construct a predictive model by integrating multiple mathematical approaches. Specifically, the Weibull model was employed to capture the non-linear dynamics of bacterial inactivation over time, while the Hill function was applied to reproduce the steep, switch-like response of bacterial survival to nisin concentration. Through this framework, we established a means to visualize and predict antimicrobial performance under diverse storage conditions in three dimensions, moving beyond static evaluations under single settings. Taken together, this study represents an important step toward both the rational design of food safety strategies and the optimized application of bacteriocins.

We selected Listeria monocytogenes as the primary target organism among Gram-positive bacteria responsible for spoilage and foodborne illness. This choice was made because L. monocytogenes has the ability to proliferate even under refrigerated conditions (around 4 °C), meaning that fresh produce stored at low temperatures can still reach infectious cell levels. Several studies have confirmed this risk; for example, Variability in Cold Tolerance of Food and Clinical Listeria monocytogenes reported that growth is possible at temperatures as low as −0.4 °C(1). In addition, contamination and recalls of ready-to-eat and refrigerated food products are frequently reported, and large outbreaks—including numerous hospitalizations and fatalities—have underscored the significant public health impact of this pathogen(2).

2 Weibull Time–Kill Model

We chose the Weibull model because it can describe non-linear survival behavior more flexibly than a first-order log-linear model. This point is important since many inactivation curves deviate from straight lines when plotted on a log scale.

For instance, Virto et al. showed that a Weibull distribution model provided an excellent description of the inactivation kinetics of Listeria monocytogenes and Escherichia coli under acid treatments, whereas a simple log-linear model could not capture the observed curvature(3). Similarly, in Listeria monocytogenes habituated on fresh produce, Poimenidou et al. (2016) fitted Weibull curves to acid challenge data and observed both shoulders and tails in survival curves depending on prior exposure, which again could not be captured by log‐linear decay(4). These findings consistently demonstrate that the Weibull approach is not just a convenient alternative but a statistically superior choice when survival data show shoulders or tails. In parallel, when analyzing antimicrobial efficacy across concentration gradients, researchers often observe sigmoidal dose–response relationships. Small doses produce little effect, but beyond a threshold, bacterial viability drops steeply before reaching a plateau. To capture this non-linear response, the Hill equation (sigmoid E–max model) is widely used. The Hill function relates concentration to effect via two intuitive parameters:

• \(\mathrm {EC50}\): the concentration at which half-maximal effect is achieved.
• \(h\) (Hill coefficient): which controls the steepness of the response curve.

For example, Rautenbach et al. (2006) fitted Hill curves to bactericidal activity data for multiple natural peptides (magainin 2, PGLa, gramicidin S, etc.), showing how both novel and well-established antimicrobials (including nisin analogs) could be quantitatively compared in terms of potency and slope(5). Likewise, similar work by others comparing natural and synthetic bacteriocins (including nisin and its variants) has found that fitting with Hill functions improves discrimination and prediction over simpler linear or threshold-based models(6).

The time-dependent reduction of bacterial counts was expressed using the Weibull model:

$$R(t)=\log_{10}(\frac{N_t}{N_0}) \;=\; - \left( \frac{t}{\delta(C, \mathrm {pH}, T)} \right)^{p}$$

• \(N(t)\): bacterial population at time \(t\)
• \(N_0\): the initial bacterial population
• \(δ(C,pH,T)\): the scale (resistance time), which depends on environment
• \(p\): the shape parameter

The effect of nisin concentration was described by the Hill function:

$$H(C)=\frac{E_{\max}C^{h}}{\mathrm {EC50}^{h}+C^{h}}\quad(0\le H\le 1)$$

• \(E(C)\): the relative antimicrobial effect (0-1)
• \(E_{max}\): the miximum achievable effect (typically set to 1)
• \(C\): the nisin concentration
• \(\mathrm {EC50}\): the concentration required to achieve 50% of the concentration-response curve
• \(h\): Hill coefficient, representing the steepness of the concentration-response curve

The environment-dependent scale parameter \(\delta (C,\mathrm {pH}, T)\) was defined as follows, where a smaller \(\delta\) indicates a faster onset of antimicrobial activity:

$$\boxed{\;\delta(C, \mathrm{pH}, T) \;=\exp (\; b_0 \;+\; b_{\mathrm{Hill}}\,H(C) \;+\; b_{\mathrm{pH2}}\,( \mathrm{pH}-\mathrm{pH}_{\mathrm{opt}} )^2 \;+\; b_T\,(T-T_{\mathrm{ref}}) \;)}$$

The parameter \(p\) characterizes the concavity or convexity of the survival curve. Model performance was evaluated using the root mean square error (RMSE).

Experimental data were obtained from the following sources:

• Listeria monocytogenes Under Acid and Antimicrobial Compounds Stress: Survival and Pathogenic Potential in Orange Juice(7)
• Biocontrol of Listeria monocytogenes on Fresh-Cut Produce by Treatment with Lytic Bacteriophages and a Bacteriocin(8)
• Nevertheless, She Resisted – Role of the Environment on Listeria monocytogenes Sensitivity to Nisin Treatment in a Laboratory Cheese Model(9)

4. Parameters

• \(R(t)\): reduction function of Listeria cells
• \(p\): shape parameter
• \(\delta\): time scale of bactericidal activity (unit: time)
• \(b_*\): regression coefficients for individual factors
• \(\mathrm {RMSE}\): magnitude of error between model prediction and experimental data (log scale)

Concentration Block

• \(b_{Hill}\): coefficient determining the strength of the concentration effect
• \(H(C)\): Hill function (0–1), representing the saturation of activity with respect to concentration \(C\)
• \(\mathrm {EC50}\): half-maximal effective concentration. A higher value indicates that a larger amount of nisin is required to achieve 50% efficacy
• \(h\): Hill coefficient, which defines the steepness of the dose–response relationship. Larger values correspond to switch-like behavior

pH Block

• \(b_{\mathrm {pH2}}\): coefficient for the squared deviation of pH
• \(\mathrm {pH}_{opt}\): optimal pH for maximum efficacy
• \(\mathrm {pH}\): actual pH of the system

Temperature Block

• \(b_T\): coefficient for storage or growth temperature
• \(T_{ref}\): reference temperature [°C], around which the effect of \(b_{\mathrm{T}}\) is evaluated
• \(T\): storage or experimental temperature [°C]

5. Parameter Estimation and Model Performance

For parameter estimation, we employed the scipy.optimize.curve_fit function to minimize the residuals between observed and predicted values using the method of least squares. Initial values and parameter bounds were specified to avoid convergence to local minima. When standard deviations were available, we applied weighted least squares, so that data points with larger experimental uncertainty were assigned smaller weights, while those with smaller error were emphasized in the fitting process.

To evaluate model performance, we calculated the root mean square error (RMSE), which reflects the magnitude of deviation between model predictions and experimental data, according to the following equation:

$$\mathrm {RMSE} = \sqrt{\frac 1 n {\sum (y_i−\hat{y}_i)^2}}$$

where \(y_i\) denotes the observed values and \(\hat{y}_i\) represents the predicted values.

Our model achieved an RMSE of 0.54 log, indicating that the average discrepancy from observations was approximately 0.5 log. Considering the inherent variability of experimental data in food microbiology, this level of accuracy can be regarded as practically acceptable.

A comparison of predicted and observed log reductions under the experimental conditions is presented in figure 1. Data points are distributed close to the diagonal, showing that the model explains the observations well. While in some cases the observed reduction was lower than predicted, large outliers were rare, suggesting robust model performance.

6. Effect of Nisin Concentration on Antimicrobial Activity

The influence of nisin concentration was incorporated into the time scale of inactivation, \(\delta\), as follows:

$$\quad \delta=\exp\{\,\cdots+\underbrace{b_{\mathrm{Hill}}H(C)}_{\text{concentration term}}+\cdots\}$$

The concentration effect itself was described by the Hill function:

$$H(C)=\frac{E_{max}C^{h}}{\mathrm {EC50}^{h}+C^{h}}\quad(0\le H\le 1)$$

• \(E(C)\): relative bactericidal effect (0–1)
• \(E_{max}\): maximum effect (typically set to 1)
• \(C\): nisin concentration
• \(\mathrm {EC50}\): concentration required to achieve 50% of the maximal effect
• \(h\): Hill coefficient, reflecting the steepness of the concentration–response relationship

In the present model, the estimated parameters describing nisin concentration dependence were \(b_{\mathrm{Hill}} \approx -10\), \(\mathrm {EC50} \approx 472\) µg/mL, and \(h \approx 5.0\). The large Hill coefficient (\(h=5\)) suggests that the antimicrobial action of nisin exhibits a switch-like response to concentration: little effect is observed below a threshold, while bactericidal activity increases sharply once the threshold is exceeded. The relatively high \(\mathrm {EC50}\) (472 µg/mL) indicates that substantial concentrations of nisin are required to achieve pronounced antimicrobial activity. This finding highlights the need for careful consideration of practical dosing levels and acceptable limits when applying nisin in food systems.

As shown in Figure 2, time–kill curves revealed that at low concentrations, bacterial counts remained nearly unchanged for extended periods. By contrast, at concentrations above the threshold (≥1000 µg/mL), a rapid log reduction was observed. This trend is consistent with the estimated parameters and visually demonstrates how the antimicrobial effect of nisin emerges only beyond a critical concentration range. Furthermore, gradual reductions were observed at all concentrations over time, with characteristic initial rapid decline followed by tailing behavior—features that are well captured by the Weibull model.

7 Effect of pH on Antimicrobial Activity

The influence of pH was incorporated into the time scale of inactivation, \(\delta\) as follows:

$$\quad \delta=\exp\{\,\cdots +\underbrace{b_{\mathrm{PH2}}(\mathrm {pH}-\mathrm {pH}_{\mathrm{opt}})^2}_{\text{pH term}}+\cdots\}$$

From parameter estimation, the values \(\mathrm {pH}_{\mathrm {opt}} \approx 6.1\) and \(b_{\mathrm{PH2}} \approx -4.2\) were obtained. In theory, the squared deviation term \((\mathrm {pH} - \mathrm {pH}_{\mathrm {opt}})^2\) should indicate that the antimicrobial effect is maximized near the optimal pH and declines as conditions deviate from this point. However, in the present model the coefficient was negative, meaning that \(\delta\) decreases (and thus the inactivation rate increases) as the system diverges from \(\mathrm {pH}_{\mathrm{opt}}\). This reversed behavior is likely influenced by experimental variability and correlations among other explanatory variables.

Nevertheless, when compared with previous findings, such as Khan et al. (2015), which reported an optimal pH range of 5–6, our results remain broadly consistent with the observation that acidic conditions enhance the antimicrobial effect. Thus, although the sign of the parameter suggests an inverted trend, the practical interpretation aligns with the literature(10).

As illustrated in Figure 3, time–kill curves varied considerably depending on pH. In particular, a pronounced log reduction was observed around pH 7, suggesting a stronger effect in alkaline conditions according to the model. By contrast, near the estimated \(\mathrm {pH}_{\mathrm {opt}}\) (\(\approx\) 5.0), reductions were relatively modest, reflecting the inverse influence of the fitted parameter. When compared with the actual pH values of fruits (Table 1), such as apples and oranges (pH 4.4–5.5 at the surface), partial discrepancies appear. Still, the general trend that antimicrobial effects are enhanced under acidic conditions is not contradicted.

8 Effect of Storage and Growth Temperature on Antimicrobial Activity

The influence of storage and growth temperature was incorporated into the inactivation time scale, \(\delta\), as follows:

$$\quad \delta=\exp\{\cdots\ +\underbrace {b_T\,(T-T_{\mathrm{ref}})}_{\text{temperature term}} +\cdots\}$$

The estimated temperature-related parameters were \(b_T \approx 0.88\), \(T_{\mathrm{ref}} \approx 7.0\). By the structure of the model, these values indicate that as storage or growth temperature rises above the reference temperature, \(\delta\) increases, resulting in a slower inactivation rate. In other words, under ambient conditions the antimicrobial effect of nisin tends to weaken, whereas at refrigeration temperatures, \(\delta\) decreases and faster bacterial inactivation can be expected. This trend reflects actual food storage practices and supports the view that antimicrobial agents are more effective at lower temperatures.

As shown in Figure 4, the time–kill curves demonstrated rapid log reduction at 3–5 °C, with substantial decreases in bacterial counts over time. By contrast, at 15 °C and 25 °C the effect was markedly weaker, with little reduction even after 24 hours. The curve at the reference temperature of 7 °C displayed intermediate behavior between refrigerated and ambient conditions, further supporting the plausibility of the parameter estimates. These results highlight that in the preservation of fresh produce, maintaining low-temperature conditions is a key factor in maximizing the antimicrobial efficacy of nisin.

9. Other Parameters

Among the estimated model parameters, the baseline term \(b_0 \approx 10\) determines the overall scale of the inactivation time \(\delta\), reflecting the tendency that a certain amount of time is required for bactericidal activity to take place. However, this term does not independently control the inactivation rate, but rather adjusts the overall balance in combination with other factors such as pH, temperature, and concentration.

The Weibull model shape parameter \(p\approx 0.13\) was found to be very small, suggesting a sharp initial reduction in bacterial counts followed by a pronounced tailing effect, where a fraction of cells persist for extended periods. This behavior is consistent with experimental observations in which resistant Listeria subpopulations remain detectable.

The overall model fit was indicated by \(\mathrm {RMSE} \approx 0.54\), corresponding to an average deviation of about 0.5 log between predictions and observations. Considering the variability commonly seen in microbial survival experiments, this level of predictive accuracy is regarded as sufficiently reliable.

As shown in Figure 5, the contour plot of log reduction across pH and nisin concentration clearly demonstrates that antimicrobial activity is maximized near pH 5.0. Conversely, under more alkaline conditions (pH ≥ 6.5), the bactericidal effect tends to plateau even with increasing concentration. This finding aligns with the estimated Hill parameters, particularly the relatively high \(\mathrm {EC50}\) (~472 µg/mL), which indicates that marked effects occur only above a certain concentration threshold.

Taken together, the model quantitatively demonstrates an important insight: nisin is most effective under acidic conditions (around pH 5.0), and at refrigeration temperatures, high concentrations can achieve reliable inactivation. These results provide supporting evidence for the practical application of nisin in food preservation.

10 Conclusion

In this section on the functional evaluation of nisin, we assessed its antimicrobial activity using a multivariate model dependent on concentration, pH, and temperature. Optimal parameters were estimated through least-squares minimization of residuals between observed and predicted values. These results clarify the environmental conditions under which nisin can function most effectively within our project framework.

In this project, the biosensor is designed to contain nisin, which can be applied to fresh produce to protect it from bacterial contamination. The analysis demonstrated that nisin acts most strongly under acidic conditions (pH \(\approx\) 5.0), is more effective at refrigeration temperatures, but requires relatively high concentrations to achieve substantial inactivation. This provides practical guidance on how to use the biosensor effectively in households and distribution settings.

With respect to pH, the findings confirm that nisin performs well in the typically acidic environments of fruits. Regarding temperature, the estimated parameter values indicate that antimicrobial activity is more sustainable under refrigeration (around 7 °C), while efficacy decreases at ambient temperatures. Although Listeria monocytogenes is capable of proliferating even under refrigerated conditions—a major food safety concern—nisin can act preventively in such environments to reduce the risk of foodborne illness. Furthermore, the model suggests that concentrations on the order of \(10^2–10^3\) µg/mL or higher are required, along with sufficient exposure time, to achieve reliable bacterial inactivation.

11 Reference

1. Myintzaw P, Pennone V, McAuliffe O, Begley M, Callanan M. Variability in Cold Tolerance of Food and Clinical Listeria monocytogenes Isolates. Microorganisms. 2022;11(1):65. doi: https://doi.org/10.3390/microorganisms11010065 ↑
2. Centers for Disease Control and Prevention (CDC). Outbreak of Listeria Infections Linked to Ready-to-Eat Food. 2025. Available at: https://www.cdc.gov/listeria/outbreaks/ready-to-eat-foods-may-2025/index.html (accessed September 28, 2025) ↑
3. Virto R, Sanz D, Álvarez I, Condon S, Raso J. Application of the Weibull model to describe inactivation of Listeria monocytogenes and Escherichia coli by citric and lactic acid at different temperatures. Journal of the Science of Food and Agriculture. 2006;86(6):865-870. doi: https://doi.org/10.1002/jsfa.2424 ↑
4. Poimenidou S, Chatzithoma D-N, Nychas G-J, Skandamis P. Adaptive Response of Listeria monocytogenes to Heat, Salinity and Low pH, after Habituation on Cherry Tomatoes and Lettuce Leaves. PLOS ONE. 2016;11(10):e0165746. doi: https://doi.org/10.1371/journal.pone.0165746 ↑
5. Rautenbach M, Gerstner GD, Vlok NM, Kulenkampff J, Westerhoff HV. Analyses of dose-response curves to compare the antimicrobial activity of model cationic alpha-helical peptides highlights the necessity for a minimum of two activity parameters. Analytical Biochemistry. 2006;350(1):81-90. doi: https://doi.org/10.1016/j.ab.2005.11.027 ↑
6. Rautenbach M, Gerstner G, Vlok M, Kulenkampff J, Westerhoff H. Analysis of concentration response curves to describe and compare the antimicrobial activity of model cationic alpha-helical peptides. Analytical Biochemistry. 2006;350(1):81-90. doi: https://doi.org/10.1016/j.ab.2005.11.027 ↑
7. Bainotti MB, Colás-Medà P, Viñas I, Alegre I. Listeria monocytogenes Under Acid and Antimicrobial Compounds Stress: Survival and Pathogenic Potential in Orange Juice. Beverages. 2025;11(4):96. doi: https://doi.org/10.3390/beverages11040096 ↑
8. Leverentz B, Conway W, Camp M, Janisiewicz W, Abuladze T, Yang M, Saftner R, Sulakvelidze A. Biocontrol of Listeria monocytogenes on Fresh-Cut Produce by Treatment with Lytic Bacteriophages and a Bacteriocin. Applied and Environmental Microbiology. 2003;69(8):4519-4526. doi: https://doi.org/10.1128/AEM.69.8.4519-4526.2003 ↑
9. Henderson LO, Erazo Flores BJ, Skeens J, Kent D, Murphy SI, Wiedmann M, Guariglia-Oropeza V. Nevertheless, She Resisted – Role of the Environment on Listeria monocytogenes Sensitivity to Nisin Treatment in a Laboratory Cheese Model. Frontiers in Microbiology. 2020;11:635. doi: https://doi.org/10.3389/fmicb.2020.00635 ↑
10. Khan A, Vu KD, Riedl B, Lacroix M. Optimization of the antimicrobial activity of nisin, Na-EDTA and pH against gram-negative and gram-positive bacteria. LWT - Food Science and Technology. 2015;61:124-129. doi: https://doi.org/10.1016/j.lwt.2014.11.035 ↑
11. Prajapati P, Porwal C, Garg M, Singh N, Sadhu SD, Chopra R, Rao ES, Agarwal A, Saeed M, Rab SO, Mahato DK, Kumar P, Kamle M, Tripathi AD. Transforming lemon peel into a sustainable reservoir of bioactives: Exploring antioxidant and anticancer potential. Food Chemistry: X. 2025;25:102172. doi: https://doi.org/10.1016/j.fochx.2025.102172 ↑
12. Pham TN, Nguyen N, Dinh D, Kieu N, Bach L, Huynh P, Nguyen M, Tran TT. Evaluate the Chemical Composition of Peels and Juice of Seedless Lemon (Citrus latifolia) Grown in Hau Giang Province, Vietnam. IOP Conference Series: Materials Science and Engineering. 2020;991:012127. doi: https://doi.org/10.1088/1757-899X/991/1/012127 ↑
13. Montero-Calderon A, Cortes C, Zulueta A, et al. Green solvents and Ultrasound-Assisted Extraction of bioactive orange (Citrus sinensis) peel compounds. Scientific Reports. 2019;9(1):16120. doi: https://doi.org/10.1038/s41598-019-52717-1 ↑
14. Henríquez C, Speisky H, Chiffelle I, Valenzuela T, Araya M, Simpson R, Almonacid S. Development of an Ingredient Containing Apple Peel, as a Source of Polyphenols and Dietary Fiber. Journal of Food Science. 2010;75(6):H172-81. doi: https://doi.org/10.1111/j.1750-3841.2010.01700.x ↑
15. Dom Z, Mujianto L, Azhar A, Masaudin S, Samsudin R. Physicochemical properties of banana peel powder in functional food products. Food Research. 2021;5:209-215. doi: https://doi.org/10.26656/fr.2017.5(S1).037 ↑
16. United States Department of Agriculture, Food and Drug Administration (FDA), Center for Food Safety and Applied Nutrition (CFSAN). Approximate pH of foods and food products. Kansas State University Research and Extension, Food Science Institute. 2007. Available at: https://cottonwood.k-state.edu/documents/fcs-docs/pH%20chart.pdf (accessed September 28, 2025) ↑