Project

< Back

Team

< Back

Members

Attributions

Behind the scenes

Model

Precision in every projection,
excellence in every decision.

Click "Learn More" to change between sections

Methods

How We Did It

Overview

The Big Picture

Results

What We Found

Overview

Introduction

References

What?

Analyze complex systems through sub-model integration of all stages of the project.

How?

Stochastic analysis of induction system, steady-state equations for purification and analysis of dose-mediated endolysin bactericidal effects.

Why?

Understand and optimize the production and effectiveness of endolysins as a promising alternative to antibiotics.

Introduction

Biological systems are highly dynamic and complex, and as we dive deeper into their intricacies, inherent variability and random fluctuations emerge. In this context, mathematical models offer a valuable tool for rigorous analysis and play a crucial role in research.¹ Thus, effectively bridging the gap between theory and practice to unlock the true potential of understanding these complex systems. Mathematical models use observation and manipulation in the same way as experiments, bypassing challenges such as high costs, time consumption, and risks.

To systematically integrate all the stages of AureoBos, from recombinant protein production to practical application, our model adopts an end-to-end Strategy that encompasses the entire process, aiming to simulate and analyze relevant stages, components, and interactions involved. viewpoint. In our approach, two sequential models are used, the first represents the upstream Describes the initial stages of the process. (induction system) and downstream Models the subsequent stages and the final outcome. (quasi-quantification process) together, while the second accounts for the implementation In vitro endolysin performance. (endolysin evaluation over target bacteria).

Initially, we use a first principle First principle models are built on a fundamental understanding of underlying ab initio physio-chemical phenomena such as mass transfer, heat transfer and mass flow. model to capture the mechanism’s behavior responsible for inducing the production of our target endolysins, we utilized isopropyl β-D-1-thiogalactopyranoside (IPTG) as an inducer, triggering the protein expression system. However, this mechanism raised important questions:

What are the specific conditions that yield the best quantitative and qualitative results in our production process?

How does the system behave under various conditions, and what factors affect its efficiency?

How robust is the system when changing the above-mentioned conditions?

The model's construction and evaluation rely on a fusion of empirical data, existing literature, and mechanistic insights. By considering multiple parameters, such as dynamic interactions among system components (kinetic coefficients/interactions), system heterogeneities, and uncertainties. We employed stochastic analysis to simulate the system's responses under various scenarios, thereby establishing a range of conditions where protein expression is most favorable. Experimental data, such as IPTG concentration and endolysin production rates, contributed to parameter estimation and model calibration.

Our model consists of a system representation through 13 ordinary differential equations (ODEs), enabling us to predict protein expression. To address the stochastic nature of the system, where behavior is influenced by probabilities² and achieve robust/realistic results, we account the uncertainty sources into our model by a Monte Carlo simulation to identify the model's parameters range and assess the sensitivity to changes. The quasi-quantification of specific endolysin yield was achieved through a grey-box Flexible and pragmatic method that acknowledges the practical limitations of accessing full internal system details while utilizing the available knowledge to create effective and efficient solutions. surrogate model, obtained combining a set of steady-state differential equations (IPTG induction mechanism), low-cost/effort lab measurements (protein distribution through SDS–PAGE) and parameter estimations of kinetic constants (Genetic Algorithm).

As mentioned, we present a combined SDS–PAGE gel electrophoresis with banding pattern recognition to enable rapid protein distribution examination. The acceleration in detection not only significantly reduces experimental costs but is also pivotal as biological system studies increasingly embrace quantitative methods, emphasizing the growing role of mathematical analysis.^3,⁴ This phase of our methodology aims to identify potential interventions or enhancements in a forward-engineering framework, advocating for the integration of a mathematical model at the conceptual design stage^5,⁶ and serves as a foundation for the (re-)design, construction, characterization, and testing of the envisioned system.

Nevertheless, a key challenge lies in accurately calibrating these models, addressing the parameter estimation problem associated with unknown rate constants for each interaction in a network.⁷ Here, allowed to determine the unknown parameters, adjusting the model to best align with experimental data, even for parameters not directly measurable, thus enhancing the model's overall fit with experimental results. Parameter robustness is crucial in modeling, as even minor variations can significantly impact the model's behavior and performance.

Moreover, having taken these crucial aspects into consideration, our focus expanded towards considering the subsequent steps required to develop a viable alternative to antibiotic therapy. This included the use of steady-state equations for purification and quantifying the endolysin concentrations, thus, the pivotal phase of our study unfolded:

How effective is the endolysin in killing the target bacteria?

Is there a dose-response relationship between endolysin concentration and bacterial kill rate?

Are there synergistic or antagonistic interactions between the endolysins?

In the culmination of our study, we sought to understand how the endolysin's mechanism interacts with various bacterial strains. This approach helped us determine if the efficacy of the endolysin is influenced by synergic or antagonistic interactions. For the analysis, we utilized a response surface Technique used to analyze the relationship between input variables and a response variable, allowing for an optimization of the response by adjusting the inputs. method (RSM) as a statistical tool in our product development work. RSM led to the identification of factor settings that collectively optimized multiple variables that align with desired responses and specifications. In our quest for optimization within our bactericidal system, RSM enabled us to make informed trade-offs between responses.

References

(1) Banwarth-Kuhn, M.; Sindi, S. How and Why to Build a Mathematical Model: A Case Study Using Prion Aggregation. J. Biol. Chem. 2020, 295 (15), 5022–5035. https://doi.org/10.1074/jbc.REV119.009851.

(2) Karlebach, G.; Shamir, R. Modelling and Analysis of Gene Regulatory Networks. Nat. Rev. Mol. Cell Biol. 2008, 9 (10), 770–780. https://doi.org/10.1038/nrm2503.

(3) O’Neill, M. A.; Denos, M.; Reed, D. Using SDS–PAGE Gel Fingerprinting to Identify Soft‐bodied Wood‐boring Insect Larvae to Species. Pest Manag. Sci. 2018, 74 (3), 705–714. https://doi.org/10.1002/ps.4766.

(4) Mathematical Modelling and Machine Learning Methods for Bioinformatics and Data Science Applications; Bianchini, M., Samponi, M. L., Eds.; MDPI: Switzerland, 2022. https://doi.org/10.3390/books978-3-0365-2841-0.

(5) Gábor, A.; Banga, J. R. Robust and Efficient Parameter Estimation in Dynamic Models of Biological Systems. BMC Syst. Biol. 2015, 9 (1), 74. https://doi.org/10.1186/s12918-015-0219-2.

(6) Arpino, J. A. J.; Hancock, E. J.; Anderson, J.; Barahona, M.; Stan, G.-B. V.; Papachristodoulou, A.; Polizzi, K. Tuning the Dials of Synthetic Biology. Microbiology 2013, 159 (Pt_7), 1236–1253. https://doi.org/10.1099/mic.0.067975-0.

(7) Mitra, E. D.; Hlavacek, W. S. Parameter Estimation and Uncertainty Quantification for Systems Biology Models. Curr. Opin. Syst. Biol. 2019, 18, 9–18. https://doi.org/10.1016/j.coisb.2019.10.006.

Methods

Underlying theory

Model building

Monte Carlo Simulation

Parameter Estimation

Response Surface Model

References

First model

Underlying theory

The IPTG induction mechanism has 3 key elements: The operator (LacO), lac repressor, and inductor (IPTG). These elements interact in a coordinated manner, as illustrated in figure 1. The process unfolds as follows:

Go over the numbers on this diagram to find more information and view the steps!

1Transcription of the LacI gene.

2Translation of the mRNA (formation of the repressor monomer).

3Formation of the repressor dimer.

4Formation of the repressor tetramer.

5The dimer binds to the LacO site (dimer-LacO complex).

6The tetramer binds to the LacO site (tetramer-LacO complex).

7 Degradation of the mRNA (repressor).

8Degradation of the repressor monomer.

9Degradation of the repressor dimer.

10Degradation of the repressor tetramer.

11Input of IPTG.

12Dissociation of dimer-LacO complex and formation of dimer-IPTG complex.

13Dissociation of tetramer-LacO complex and formation of tetramer-IPTG complex.

14Formation of dimer-IPTG complex.

15Formation of tetramer-IPTG complex.

16Degradation of the dimer-IPTG complex.

17Degradation of the tetramer-IPTG complex.

18Transcription of the endolysin gene.

19Translation of the mRNA (formation of endolysin).

20Degradation of the mRNA (endolysin).

21Degradation of the endolysin.

Figure 1. IPTG induction steps when Escherichia coli is found under conditions with IPTG and no glucose.

Production of Lac repressor. Steps 1-4

The system's dynamics are under the influence of the Lac repressor, which is a regulatory molecule encoded by the LacI regulatory gene (G). This protein exhibits constitutive expression, forming monomers (R_M), dimers (R_D), and tetramers (R_T) through subunit bonding.³ The chemical reaction for modeling the repressor mRNA production is assumed to be governed by the sole dependence on the quantity of the LacI gene. Additionally, the production of monomer is dependent on the mRNA of repressor (M_LI) concentration. The favored formation is Tetramers over dimers over monomers.

G \overset{k_{1}}{\to} M_{L I}

(i)

M_{L I} \overset{k_{2}}{\to} M_{L I} + R_{M}

(ii)

2 R_{M} ⇄_{k_{4}}^{k_{3}} R_{D}

(iii,iv)

2 R_{D} ⇄_{k_{6}}^{k_{5}} R_{T}

(v,vi)

Formation of repressor-Operator complex. Steps 5-6

In its tetrameric and dimeric configurations, the repressor exhibits strong affinity for the operator (O). This binding event effectively obstructs RNA polymerase’s access to the promoter region, resulting in the inhibition of endolysin gene transcription. Consequently, in the absence of other factors, endolysin production is low or absent (leak transcript). On the contrary, the monomeric form of the lac repressor demonstrates a reduced affinity for the operator and lacks the capability to effectively bind to it.⁴

R_{D} + O ⇄_{k_{8}}^{k_{7}} R_{D} O

(vii,viii)

R_{T} + O ⇄_{k_{10}}^{k_{9}} R_{D} O

(ix,x)

Input of IPTG. Step 11

The entry of IPTG into E. coli BL21 cells occurs through two potential mechanisms: simple diffusion and active transport facilitated by the lactose permease, encoded by the LacY gene. The active transport pathway has been reported to not be significant on E. coli, thus, we won't take it into account in our model.⁵ As a result, we focus solely on the predominant method of IPTG entry, which is simple diffusion. In this process, IPTG molecules passively move across the cell membrane in response to the concentration gradient, without the need for energy expenditure or assistance from specific transporters.

We employ Fick's law of diffusion (Eq. 1) to characterize the process of simple diffusion. The equation takes into account 2 variables: IPTG concentration and absorption surface. Notably, in our model, we maintain the absorption surface as a constant, disregarding any membrane expansion during binary fission. This choice enables us to establish all our mathematical equations as dependent solely on IPTG concentration. In this context, where Δx denotes the membrane thickness, 'K' represents the diffusion constant, 'A' stands for the absorption surface, 'IPTG_E ' signifies the IPTG concentration in the culture medium, and 'IPTG_I' refers to the IPTG concentration within the bacterial cells.

\frac{Δ I P T G_{I}}{Δ t} = K A \frac{Δ I P T G_{E}}{Δ x}

(1)

Formation of repressor-IPTG complex. Steps 12-15

Upon the introduction of IPTG into the system, it binds with the dimeric and tetrameric forms of the repressor protein. It’s important to note that each monomer within both the dimer and tetramer structure can accommodate a single IPTG molecule for binding. Upon binding, IPTG triggers an allosteric structural transformation in the tetrameric repressor, resulting to a decrease in its affinity for the operator.⁶ Consequently, the monomeric repressor loses its capacity to effectively bind to the operator, allowing RNA polymerase’s access to the promoter. This, in turn, enables the transcription and subsequent production of the endolysins.

The formation of the IPTG-repressor complex can occur through two mechanisms:³

The repressor exists in a free state.
The repressor is part of the repressor-operator complex, where it is already bound to the operator sequence.

When the repressor is bound to the operator, the introduction of IPTG interferes with the repressor-operator complex, leading to its disruption. As a result, the extent of repressor inactivation is directly correlated with the number of binding sites occupied by IPTG. Interestingly, even with just two IPTG molecules bound to specific sites, the repressor can undergo effective inactivation.⁷ Figure 2 illustrates the different conformations that can exist.

We are considering these distinct conformations to deactivate the tetramer repressor, which requires a Hill coefficient of 2-3. However, it's noteworthy that only two IPTG molecules are required to inactivate the dimeric repressor.

2 I + R_{D} ⇄_{k_{12}}^{k_{11}} I_{2} R_{D}

(xi,xii)

2 I + R_{D} O ⇄_{k_{14}}^{k_{13}} I_{2} R_{D} + O

(xiii,xiv)

n I + R_{T} ⇄_{k_{16}}^{k_{15}} I_{n} R_{T}

(xv,xvi)

n I + R_{T} O ⇄_{k_{18}}^{k_{17}} I_{n} R_{T} + O

(xvii,xviii)

Protein production. Steps 18-19

As the operator becomes unoccupied due to the inactivation of the lac repressor through IPTG binding, RNA polymerase gains the ability to initiate transcription of the endolysin gene. In our model, the transcription of endolysin is solely contingent on the concentration of free operator sites.

E. coli recombinant expression systems that utilize lac operon control elements to modulate gene expression, low levels expression are produced during bacterial growth prior to induction, known as leakage. This leaking limits cell growth when toxic proteins for E. coli are produced and cause instability in the pET expression system.⁸

The concentration of recombinant proteins in bacterial cells is crucial for large-scale production, as high levels can be toxic, reducing yields. However, for the sake of model simplicity and tractability, we have omitted toxicity considerations.

O \overset{k_{19}}{\to} O + M_{E}

(xix)

R_{D} O \overset{k_{20}}{\to} R_{D} O + M_{E}

(xx)

R_{T} O \overset{k_{21}}{\to} R_{T} O + M_{E}

(xxi)

M_{E} \overset{k_{22}}{\to} M_{E} + E

(xxii)

Degradation of mRNA and protein compounds. Steps 7-10, 16-17, 20-21

The non-hydrolyzable chemical bond in IPTG, specifically involving the sulfur (S) atom, is a critical aspect of its mechanism of action within the cell. This unique characteristic prevents IPTG from being metabolized or degraded by cellular machinery. In our model, we take into account this specific feature, making an important assumption that the compounds that form complexes with IPTG are subject to degradation, leaving behind the free IPTG molecules that were involved in the complex formation.

It is worth noting that when the lac repressor exists in its dimeric or tetrameric form and is bound to the operator, it is rendered non-degradable. This is because these bound repressor complexes are part of the regulatory mechanism that inhibits the transcription of the endolysin gene. As such, the bound dimers and tetramers are not subject to degradation in the model.

M_{L I} \overset{k_{23}}{\to} \emptyset

(xxiii)

R_{M} \overset{k_{24}}{\to} \emptyset

(xxiv)

R_{D} \overset{k_{25}}{\to} \emptyset

(xxv)

R_{T} \overset{k_{26}}{\to} \emptyset

(xxvi)

I_{2} R_{D} \overset{k_{27}}{\to} \emptyset

(xxvii)

I_{n} R_{T} \overset{k_{28}}{\to} \emptyset

(xxviii)

M_{E} \overset{k_{29}}{\to} \emptyset

(xxix)

E \overset{k_{30}}{\to} \emptyset

(xxx)

Model building

Chemical reactions

Drawing from the available literature, we have acquired the essential chemical reactions governing our system's behavior. These reactions laid the groundwork for formulating our ODEs.

Table 1. Chemical reactions on the IPTG induction model.

Number	Chemical reactions	Meaning
(1)	N/A	Entrance of IPTG into E. coli (Simple diffusion)
(i)	$M_{L I} \overset{k_{2}}{\to} M_{L I} + R_{M}$	Production of the repressor mRNA.
(ii)	$M_{L I} \overset{k_{2}}{\to} M_{L I} + R_{M}$	Translation of the repressor mRNA (monomer formation)
(iii,iv)	$2 R_{M} ⇄_{k_{4}}^{k_{3}} R_{D}$	Formation of dimers (repressor).
(v,vi)	$2 R_{D} ⇄_{k_{6}}^{k_{5}} R_{T}$	Formation of tetramers (repressor).
(vii,viii)	$R_{D} + O ⇄_{k_{8}}^{k_{7}} R_{D} O$	Binding of the dimer to LacO (operator).
(ix,x)	$R_{T} + O ⇄_{k_{10}}^{k_{9}} R_{T} O$	Binding of the tetramer to operator.
(xi,xii)	$2 I + R_{D} ⇄_{k_{12}}^{k_{11}} I_{2} R_{D}$	Binding of IPTG to the dimer.
(xiii,xiv)	$2 I + R_{D} O ⇄_{k_{14}}^{k_{13}} I_{2} R_{D} + O$	Binding of IPTG to the dimer-operator complex.
(xv,xvi)	$n I + R_{T} ⇄_{k_{16}}^{k_{15}} I_{n} R_{T}$	Binding of IPTG to tetramer.
(xvii,xviii)	$n I + R_{T} O ⇄_{k_{18}}^{k_{17}} I_{n} R_{T} + O$	Binding of IPTG to the tetramer-operator complex.
(xix)	$O \overset{k_{19}}{\to} O + M_{E}$	Transcription of endolysin.
(xx)	$R_{D} O \overset{k_{20}}{\to} R_{D} O + M_{E}$	Transcript leak from endolysin in dimers.
(xxi)	$R_{T} O \overset{k_{21}}{\to} R_{T} O + M_{E}$	Transcript leak from endolysin in tetramers.
(xxii)	$M_{E} \overset{k_{22}}{\to} M_{E} + E$	Translation of the endolysin.
(xxiii)	$M_{L I} \overset{k_{23}}{\to} \emptyset$	Degradation of the repressor mRNA.
(xxiv)	$R_{M} \overset{k_{24}}{\to} \emptyset$	Degradation of the repressor monomers.
(xxv)	$R_{D} \overset{k_{25}}{\to} \emptyset$	Degradation of the repressor dimers.
(xxvi)	$R_{T} \overset{k_{26}}{\to} \emptyset$	Degradation of the repressor tetramers
(xxvii)	$I_{2} R_{D} \overset{k_{27}}{\to} \emptyset$	Degradation of the dimer-IPTG complex.
(xxviii)	$I_{n} R_{T} \overset{k_{28}}{\to} \emptyset$	Degradation of the tetramer-IPTG complex.
(xxix)	$M_{E} \overset{k_{29}}{\to} \emptyset$	Degradation of the endolysin mRNA.
(xxx)	$E \overset{k_{30}}{\to} \emptyset$	Degradation of the endolysin.

Ordinary Differential Equations

The biochemical reactions in the system were captured through the development of an ordinary differential equations (ODEs) model. All reactions within this model were defined using the law of mass action, offering a dynamic and quantitative description of the system. The ODEs express how the system's temporal concentrations change with time, providing a framework to analyze and predict its behavior.

The model was implemented in MATLAB with the ODE15s solver.

ODEs

$\frac{d [[M_{L I}]]}{d t} = k_{2} [[G, ]] - k_{24} [[M_{L I}]]$

(2)

\frac{d [[R_{M}]]}{d t} = k_{3} [[M_{L I}]] - 2 ((, k_{4} {[[R_{M}]]}^{2} - k_{5} [[R_{D}])]) - k_{25} [[R_{M}]]

(3)

\frac{d [[R_{D}]]}{d t} = ((, k_{4} {[[R_{M}]]}^{2} - k_{5} [[R_{D}])]) - 2 ((, k_{6} {[[R_{D}]]}^{2} - k_{7} [[R_{T}, ], )]) - ((k_{8} [[R_{D}]] [[, O, ]] - k_{9} [[R_{D} O], ], )) -

((, k_{12} [{[I, ]]}^{2} [[R_{D}, ]] - k_{13} [[I_{2} R_{D}, ]] - k_{26} [[R_{D}, ], )])

(4)

\frac{d [[R_{T}]]}{d t} = ((k_{6} {[[R_{D}]]}^{2} - k_{7} [[R_{T}, ], )]) - ((, k_{10} [[R_{T}], [] [O, ]] - k_{11} [[R_{T} O], ], )) - ((, k_{16} {[[, I, ]]}^{n} [[R_{T}]] - k_{17} [[I_{n} R_{T}, ], )]) - k_{27} [[R_{T}]]

(5)

\frac{d [[, O, ]]}{d t} = ((k_{8} [[R_{D}]] [[, O, ]] - k_{9} [[R_{D} O], ], )) - ((k_{10} [[R_{T}]] [[, O, ]] - k_{11} [[R_{T} O], ], )) + ((k_{14} {[[, I, ]]}^{2} [[R_{D} O, ]] - k_{15} [[I_{2} R_{D}]] [[, O, ], )]) +

((k_{18} {[[, I, ]]}^{n} [[R_{T} O, ]] - k_{19} [[I_{n} R_{T}]] [[, O, ], )])

(6)

\frac{d [[R_{D} O, ]]}{d t} = ((k_{8} [[R_{D}]] [[, O, ]] - k_{9} [[R_{D} O], ], )) - ((k_{14} {[[, I, ]]}^{2} [[R_{D} O, ]] - k_{15} [[I_{2} R_{D}]] [[, O, ], )])

(7)

\frac{d [[R_{T} O, ]]}{d t} = ((, k_{10} [[R_{T}], [] [O, ]] - k_{11} [[R_{T} O], ], )) - ((k_{18} {[[, I, ]]}^{n} [[R_{T} O, ]] - k_{19} [[I_{n} R_{T}]] [[, O, ], )])

(8)

\frac{d [[, I, ]]}{d t} = - 2 ((, k_{12} [{[I, ]]}^{2} [[R_{D}, ]] - k_{13} [[I_{2} R_{D}, ], )] - 2 ((k_{14} {[[, I, ]]}^{2} [[R_{D} O, ]] - k_{15} [[I_{2} R_{D}]] [[, O, ], )]) - n ((, k_{16} {[[, I, ]]}^{n} [[R_{T}]] - k_{17} [[I_{n} R_{T}, ], )]))

- n ((k_{18} {[[, I, ]]}^{n} [[R_{T} O, ]] - k_{19} [[I_{n} R_{T}]] [[, O, ], )]) + 2 k_{28} [[I_{2} R_{D}]] + 2 k_{29} [[I_{n} R_{T}]]

(9)

\frac{d [[I_{2} R_{D}]]}{d t} = ((, k_{12} [{[I, ]]}^{2} [[R_{D}, ]] - k_{13} [[I_{2} R_{D}, ], )] + ((k_{14} {[[, I, ]]}^{2} [[R_{D} O, ]] - k_{15} [[I_{2} R_{D}]] [[, O, ], )]) - k_{28} [[I_{2} R_{D}]])

(10)

\frac{d [[I_{n} R_{T}]]}{d t} = ((, k_{16} {[[, I, ]]}^{n} [[R_{T}]] - k_{17} [[I_{n} R_{T}, ], )]) + ((k_{18} {[[, I, ]]}^{n} [[R_{T} O, ]] - k_{19} [[I_{n} R_{T}]] [[, O, ], )]) - k_{29} [[I_{n} R_{T}]]

(11)

\frac{d [[M_{E}]]}{d t} = k_{20} [[, O, ]] + k_{21} [[R_{D} O, ]] + k_{22} [[R_{T} O]] - k_{30} [[M_{E}]]

(12)

\frac{d [[, E, ]]}{d t} = k_{23} [[M_{E}]] - k_{31} [[, E, ]]

(13)

Table 2. Constants on the IPTG induction model

Constant	Meaning	Value		Unit	Reference
Constant	Meaning	LysK	LysCSA13	Unit	Reference
$k_{1}$	Simple diffusion constant of IPTG*	0.92		$m i n^{- 1}$	⁴
$k_{2}$	Transcription rate of LacI*	0.23		$n M \cdot m i n^{- 1}$	⁴
$k_{3}$	Translation rate of LacI*	15		$m i n^{- 1}$	⁴
$k_{4}$	Dimerization rate of LacI*	50		$n M \cdot m i n^{- 1}$	⁴
$k_{5}$	LacI dimerization rate constant*	10^-3		$m i n^{- 1}$	⁴
$k_{6}$	Tetramerization rate of LacI	172806.1332	528929.5329	$n M^{- 1} \cdot m i n^{- 1}$	estimated
$k_{7}$	Dissociation rate of LacI tetramers	0.004646925019	0.004682362659	$m i n^{- 1}$	estimated
$k_{8}$	Association rate of the first repression mechanism*	960		$n M^{- 1} \cdot m i n^{- 1}$	⁴
$k_{9}$	Dissociation rate of the first repression mechanism*	2.4		$m i n^{- 1}$	⁴
$k_{10}$	Association rate of the second repression mechanism	8289.871938	5870.331195	$n M^{- 1} \cdot m i n^{- 1}$	estimated
$k_{11}$	Dissociation rate of the second repression mechanism	1.401669318	1.596701205	$m i n^{- 1}$	estimated
$k_{12}$	Association rate of the second derepression mechanism*	3x10^-7		$n M^{- 2} \cdot m i n^{- 1}$	⁴
$k_{13}$	Dissociation rate of the first derepression mechanism*	12		$m i n^{- 1}$	⁴
$k_{14}$	Association rate of the second derepression mechanism*	3x10^-7		$n M^{- 2} \cdot m i n^{- 1}$	⁴
$k_{15}$	Dissociation rate of the second derepression mechanism*	4.8x10^-7		$n M^{- 1} \cdot m i n^{- 1}$	⁴
$k_{16}$	Association rate of the third derepression mechanism	7.36E-07	3.94E-07	$n M^{- 2} \cdot m i n^{- 1}$	estimated
$k_{17}$	Dissociation rate of the third derepression mechanism	2988.006153	2284.487239	$m i n^{- 1}$	estimated
$k_{18}$	Association rate of the fourth derepression mechanism	9.87E-08	1.00E-07	$n M^{- 2} \cdot m i n^{- 1}$	estimated
$k_{19}$	Dissociation rate of the fourth derepression mechanism	2411.253642	1960.781346	$m i n^{- 1}$	estimated
$k_{20}$	Transcriptional rate of endolysin	0.1361368059	0.1283765476	$m i n^{- 1}$	estimated
$k_{21}$	Transcriptional leak rate of endolysin in dimers	0.07776261945	0.2261425676	$m i n^{- 1}$	estimated
$k_{22}$	Transcriptional leak rate of endolysin in tetramers	0.4241461762	0.3737824434	$m i n^{- 1}$	estimated
$k_{23}$	Translation rate of the endolysin	13.7090742	10.93457685	$m i n^{- 1}$	estimated
$k_{24}$	Degradation rate constant of the repressor mRNA*	0.462		$m i n^{- 1}$	⁴
$k_{25}$	Degradation rate constant of the repressor monomer*	0.2		$m i n^{- 1}$
$k_{26}$	Degradation rate constant of the repressor dimer*	0.2		$m i n^{- 1}$	⁴
$k_{27}$	Degradation rate constant of the repressor tetramer	0.08567815958	0.1271362511	$m i n^{- 1}$	estimated
$k_{28}$	Degradation rate constant of the dimer-IPTG complex*	0.2		$m i n^{- 1}$	⁴
$k_{29}$	Degradation rate constant of the tetramer-IPTG complex	0.05925751146	0.1315557746	$m i n^{- 1}$	estimated
$k_{30}$	Degradation rate constant of the repressor monomer	0.05214786934	0.1212103757	$m i n^{- 1}$	estimated
$k_{31}$	Degradation rate constant of endolysin mRNA	0.8233674927	1.531937249	$m i n^{- 1}$	estimated

Table 3. Initial conditions on the IPTG induction model.

Variable	Meaning	Value	Unit
$I$	Free IPTG inside the cell	0	$n M$
$I_{e x t}$	Extracellular IPTG	100	$n M$
$O$	Free operator	50	$n M$
$G$	Gene of the repressor	10	$n M$
$M_{L I}$	mRNA of the repressor	50	$n M$
$R_{M}$	Monomer of the repressor	50	$n M$
$R_{D}$	Dimer of the repressor	10	$n M$
$R_{T}$	Tetramer of the repressor	40	$n M$
$R_{D} O$	Dimer-operator complex	15	$n M$
$I_{2} R_{D}$	Dimer-IPTG complex	0	$n M$
$R_{T} O$	Tetramer-operator complex	45	$n M$
$I_{n} R_{T}$	Tetramer-IPTG complex	0	$n M$
$M_{E}$	Endolysin mRNA	5	$n M$
$E$	Endolysin	3	$n M$
$\emptyset$	Degradation	N/A	$n M$

Monte Carlo Simulation

Precise predictions of biological processes rely on understanding intricate networks. Nonetheless, this understanting is hindered by the lack of data concerning numerous model parameters, especially reaction kinetics. In our system of ODEs, we encountered a shorfall of associated parameters for some of the reaction kinetics presented in Table 2. Selecting the suitable -range of- values presented a challenge as it should encompass a broad spectrum to ensure the experimental results are representative, detailed enough to capture essential information, but not overly extensive that it becomes too complex and time-consuming due to limitations in computational resources.

Monte Carlo (MC) simulations have been employed to predict the effects of stochasticity among reaction events, as well as in model developments where parameters are unknown, and for parameter estimation.^15-18 Furthermore, MC simulations demonstrate a noteworthy level of accuracy, efficiency in computation, and ease of implementation. In our model, we employ a MC simulation to conduct 10,000 iterative excecutions utilizing randomly assing values for the unknown kinetic constants in Table 2.

Through this approach, a diverse and extensive array of potential parameter sets is generated, facilitating a probabilistic exploration of the parameter space. Subsequently, we leverage the established range of viable kinetic constant values obtained from the MC simulation as an input for a machine-learning-based parameter estimation algorithm, aiming to narrow down said parameter range.

Parameter estimation

Subsequently, employing the collected data set of sampled parameters, we proceed to estimate the specific parameters tied to the kinetic aspects of the system. We present a computational approach for the estimation of kinetic constants, focusing on the induction phase of protein expression, through the application of a Genetic Algorithm (GA).

Combining the range of parameters obtained through MC simulation with a GA -or any suitable parameter estimation algorithm- is essential when modeling systems characterized by indeterminate kinetic constants. The empirical data sourced from SDS-PAGE gel electrophoresis was analyzed through ImageJ. The initial preprocessing phase involves the segmentation of gel lanes and accurate delineation of protein bands, extracting quasi-quantitative concentration data across different time points. This processed concentration data forms the foundation for constructing an objective function within the GA, to iteratively optimize the kinetic parameters.²⁰ The objective function evaluates the fitness of each parameter set by comparing simulated protein concentration profiles with the experimentally derived data. Employing a generational approach, the GA refines the kinetic parameter space, converging towards optimal or near-optimal solutions that exhibit the highest fidelity with observed protein distribution.

Second model

Response Surface Method

Background

To address the preparation of an intramammary infusion composed of our three endolysins LysK-ABD-SH3B30, LysCSA13-ABD, and PCNP-CecA-LysSS (which will further be referred as LysK, LysCSA13 and LysSS), it was imperative to determine the optimal composition (combination of concentrations) for targeting each specific bacterium. Initially, experiments were envisioned utilizing purified recombinant proteins; however, challenges were encountered in the purification process. Consequently, a decision was made to conduct bactericidal assays using crude extracts of each endolysin. Despite this modification, the primary objective remains to attain the most effective formulation of AureoBos.

The presence of endolysins targeting different bonds within the peptidoglycan suggests a promising potential for a synergistic effect against their target bacteria when these enzymes are used in combination.²⁰ While one-factor-at-a-time approach is commonly used to optimize components and conditions, it is time-consuming and tends to overlook the interactions between independent variables.²¹ To screen and analyze interactions between the endolysins, particularly focusing on their combined effect, we aim to analyze the potential synergistic impact, finding the optimum quickly and efficiently.

Response Surface Methodology (RSM) is a surrogate modelling technique that combines mathematical and statistical techniques valuable for construction models and analyzing scenarios where multiple independent variables or controllable factors influence a dependent variable or response. RSM efficiently extracts a wealth of information through minimal effort.²² As such, we utilize RSM to estimate a response based on a combination of factor levels and determining the optimum concentratios to maximize the endolysins’ performance. The optimum concentrations were first determined by screening independent factors (protein concentrations of each construct) and selecting desired responses (bactericidal effect), followed by choosing the experimental design strategy and running the experiments to obtain the results. After that, the model is confirmed by utilizing response graphs and analysis of variance, and optimal conditions are determined in the end.

Experimental set-up

In our case, the three factors, namely LysK, LysCSA13, and LysSS, in three levels within 875-2625 µg mL^-1, 500-1500 µg mL^-1, 25-75 µg mL^-1 respectively, were assessed using three-level full factorial design. Design that involves systematically varying all possible combinations of factor levels to study the impact of multiple factors and their interactions on the response. The concentration levels were determined based on total protein concentration in each crude extract, categorizing them as low, middle, and high boundaries. Our model operates under the assumption that the optimal conditions fall within the specified ranges. The focus was solely on assessing the efficacy of the endolysins, thus, we maintained constant temperature and pH conditions throughout the experiments. Each experiment was conducted in triplicate.

Table 4. Values of protein concentrations in µg/mL for each variable.

		Level
Variable	-1	0	1
LysK	875	1750	2625
LysCSA13	500	1000	1500
LysSS	25	50	75

As previously mentioned, we examined combinations of varying levels of protein concentration in each crude extract. However, due to constraints, these combinations were tested against only three of the four target bacteria (E. coli, S. aureus, S. agalactie). The bactericidal assay involved the addition of the combination crude extracts of LysK, LysCSA13, and LysSS to the respective bacteria. In all combinations, we employed resazurin to evaluate cell viability post-treatment with crude extracts, measuring fluorescence at 560 nm excitation and 590 nm emission. The full description of our methods can be found in the Experiments page. Additionally, we established a correlation between the measured fluorescence and the population of viable cells, expressed through the following equation:

B a c t e r i c i d a l a c t i v i t y % = 100 - (\frac{[s a m p l e]}{P o s i t i v e C o n t r o l - N e g a t i v e C o n t r o l} \cdot l 00)

Where cell viability (% control) represents the proportion of viable cells remaining post-treatment. A sample refers to a specific concentration combination being assayed. The positive control is the corresponding bacteria treated with resazurin, while the negative control comprises the crude extract of E. coli BL21(DE3), the corresponding bacteria and resarzurin.

Data analysis

The data collected from our experiments (available in the datatables 1 ,2,3) was analyzed using the Design Expert software by Stat-Ease. Our objective was to minimize the number of viable cells of each bacteria, and for this purpose, we utilized the software to input the calculated correlations for each combination of the tested factors. The following methodology was applied to each case:

Creating appropriate model

In order to construct an accurate model from our data, it's crucial to have a normally distributed dataset. We achieved this by implementing a Box-Cox Power transformation to find the best λ that maximizes the log-likelihood of the transformed data, making it as close to normally distributed as possible. plot. After applying the suitable transformation, we confirmed the normal distribution through a normal probability plot and determined the best-fit equation recommended by Design Expert out of three options: linear, 2FI, or quadratic. If the p-value exceeded 0.05 and the term wasn't necessary hierarchically, it was eliminated. Below is a summary outlining the models for each bacterium.

S. aureus

For the S. aureus response, the cuadratic model was recommended. Box-Cox plot for the crude extracts data for bactericiday activity (%) against S. aureus. Based on this, the software recommended not performing a transformation.

Figure 3.Box-Cox plot for S. aureus data. The green line represents the best λ (2.77). The red line represent the confidence interval (-0.07, 5.38), and the blue line represents the current transformation (1).

Based on Figure 4, it is evident that the data closely follows a normal distribution. The points on the plot align well with the straight line, indicating that the distribution of the data is approximately normal. This conformity to a normal distribution is important for various statistical analyses that assume normality, as it ensures the validity and reliability of the results derived from those analyses. It signifies that the underlying assumptions of these statistical methods, based on normality, are met by the dataset.

Figure 4. Normal probability plot of S. aureus data. No transformation applied
ANOVA, or Analysis of Variance, is a statistical method used to assess differences in means among two or more groups or treatments. It aims to determine whether these differences are significant or if they could likely be due to random chance. The key values in ANOVA provide crucial insights into the variability within and between groups. The "Sum of Squares" (SS) quantifies the dispersion of data, categorizing it into between-group and within-group variance. "Degrees of Freedom" (df) represent the number of values that are free to vary within the statistical analysis. "Mean Square" (MS) is a measure of variance, enabling the computation of the F-value—a ratio that assesses the significance of group differences. The associated p-value helps researchers conclude whether the observed differences are statistically meaningful or likely a result of chance.

Table 5. ANOVA table of the bactericidal activity agains S. aureus. Quadratic equation: no transformation applied.

Source	Sum of squeares	df	Mean square	F-value	p-value
Model	1087.83	9	120.87	5.36	0.0015	significant
A-LysK	218.41	1	2018.41	9.69	0.0063
B-LysCSA13	368.74	1	368.74	16.36	0.0008
C-LysSS	5.86	1	5.86	0.2600	0.6167
AB	203.78	1	103.78	4.61	0.0466
AC	136.22	1	136.22	6.04	0.0250
BC	2.00	1	2.00	0.0888	0.7693
A²	16.01	1	16.01	0.7103	0.4110
B²	203.35	1	203.35	9.02	0.0080
C²	33.46	1	33.46	1.49	0.2396
Residual	383.09	17	22.53
Cor total	1470.92	26

The terms BC, A², and C² having p-values larger than 0.05 indicate that these terms are not statistically significant in the model. In statistical analysis, a p-value greater than 0.05 suggests that there is insufficient evidence to reject the null hypothesis, which often means the variable or term is not contributing significantly to the model's explanatory power. Thus, we opted to remove them from the model.

Table 6. ANOVA table of the bactericidal activity agains S. aureus. Quadratic equation without the terms BC, A², and C²: no transformation applied.

Source	Sum of squeares	df	Mean square	F-value	p-value
Model	1036.36	6	172.73	7.95	0.0002	significant
A-LysK	218.41	1	218.41	10.05	0.0048
B-LysCSA13	368.74	1	368.74	16.97	0.0005
C-LysSS	5.86	1	5.86	0.2697	0.6092
AB	103.78	1	103.78	4.78	0.0409
AC	136.22	1	136.22	6.27	0.0211
B²	203.35	1	203.35	9.36	0.0062
Residual	434.56	20	21.73
Cor total	1470.92	26

S. agalactiae

For the S. agalactie response, a two-factor (2FI) model was recommended. Box-Cox plot for the crude extracts data for bactericidal activity (%) against S. agalactie. Based on this, the software recommended not performing a transformation.

Figure 5. Box-Cox plot for S. agalactie data. The green line represents the best λ (0.29). Confidence interval (-4.32, 4.90), and the blue line represents the current transformation (1).

As depicted in Figure 6, the dataset also exhibits a notable alignment of points along the plotted straight line, suggesting an approximate normal distribution of the data. Thus, a transformation was not necessary.

Figure 6. Normal probability plot of S. agalactie data. No transformation applied.

Table 7. ANOVA table of the bactericidal activity agains S. agalactie. 2FI equation: no transformation applied.

Source	Sum of squares	df	Mean square	F-value	p-value
Model	375.59	6	62.60	5.58	0.0015	significant
A-LysK	82.65	1	82.65	7.37	0.0133
B-LysCSA13	0.4640	1	0.4640	0.0414	0.8408
C-LysSS	2.18	1	2.18	0.1942	0.6642
AB	2.15	1	2.15	0.1918	0.6661
AC	286.16	1	286.16	25.53	<0001
BC	1.99	1	1.99	0.1778	0.6778
Residual	224.20	20	11.21
Cor Total	599.79	26

As indicated in Table 8, the sole term displaying statistical significance is AC. Consequently, we made the decision to refine our model, focusing solely on the AC term.

Table 8. ANOVA table of the bactericidal activity against S. agalactie. 2FI equation without AB and BC: no transformation applied.

Source	Sum of squares	df	Mean square	F-value	p-value
Model	370.99	3	123.66	12.43	<0.0001	significant
A-LysK	82.65	1	82.65	8.31	0.0084
C-LysSS	2.18	1	2.18	0.2188	0.6443
AC	286.16	1	286.16	28.77	<0.0001
Residual	228.80	23	9.95
Cor Total	599.79	26

E. coli

For the E. coli response, linear and quadratic models were recommended. Box-Cox plot for the crude extracts data for bactericiday activity (%) against E. coli. Based on this, the software recommended not performing a transformation.

Figure 7.Box-Cox plot for E. coli data. The green line represents the best λ (-0.77). The red lines represent the confidence interval (-3.11, 1.23), and the blue line represents the current transformation (1).
Illustrated in Figure 8, the dataset displays a significant alignment of data points along the plotted straight line, implying an approximate normal distribution of the data. Thus, a transformation was also deemed unnecessary.

Figure 8.Normal probability plot of E. coli data. No transformation applied.

Table 9. ANOVA table of the bactericidal activity against E. coli. Quadratic equation: no transformation applied.

Source	Sum of squares	df	Mean square	F-value	p-value
Model	685.35	9	76.15	3.35	0.0153	significant
A-LysK	3.87	1	3.87	0.1703	0.6850
B-LysCSA13	229.34	1	229.34	10.09	0.0055
C-LysSS	90.14	1	90.14	3.96	0.0628
AB	18.60	1	18.60	0.8180	0.3784
AC	6.41	1	6.41	0.2819	0.6023
BC	54.87	1	54.87	2.41	0.1387
A²	19.62	1	19.62	0.8629	0.3659
B²	195.97	1	195.97	8.62	0.0092
C²	66.53	1	66.53	2.93	0.1053
Residual	386.55	17	22.74
Cor Total	1071.90	26

As outlined in Table 10, the sole significant term appears to be B². However, this finding contradicts the expected outcome, considering our target is E. coli, a gram-negative bacterium. Given that LysCSA13 is specific to S. aureus, it would not exhibit significant bactericidal activity in this particular experiment. As such, we opted to leave the model to explore the potential effects of the C terms, even if they were initially deemed insignificant.

Table 10. ANOVA table of the bactericidal activity against E. coli. Quadratic equation without A, AB, AC, BC, and A²: no transformation applied.

Source	Sum of squares	df	Mean square	F-value	p-value
Model	581.98	4	145.49	6.53	0.0013	significant
B-LysCSA13	229.34	1	229.34	10.30	0.0040
C-LysSS	90.14	1	90.14	4.05	0.0566
B²	195.97	1	195.97	8.80	0.0071
C²	66.53	1	66.53	2.99	0.0979
Residual	489.93	22	22.27
Cor Total	1071.90	26

References

(1) Lewis, M.; Chang, G.; Horton, N. C.; Kercher, M. A.; Pace, H. C.; Schumacher, M. A.; Brennan, R. G.; Lu, P. Crystal Structure of the Lactose Operon Repressor and Its Complexes with DNA and Inducer. Science 1996, 271 (5253), 1247–1254. https://doi.org/10.1126/science.271.5253.1247.

(2) Stamatakis, M.; Mantzaris, N. V. Comparison of Deterministic and Stochastic Models of the Lac Operon Genetic Network. Biophys. J. 2009, 96 (3), 887–906. https://doi.org/10.1016/j.bpj.2008.10.028.

(3) Hansen, L. H.; Knudsen, S.; Sørensen, S. J. The Effect of the lacY Gene on the Induction of IPTG Inducible Promoters, Studied in Escherichia Coli and Pseudomonas Fluorescens. Curr. Microbiol. 1998, 36 (6), 341–347. https://doi.org/10.1007/s002849900320.

(4) Merino, J.; Noriega, M. Transporte de Membrana, 2011. https://ocw.unican.es/course/view.php?id=44#section-0.

(5) Taraban, M.; Zhan, H.; Whitten, A. E.; Langley, D. B.; Matthews, K. S.; Swint-Kruse, L.; Trewhella, J. Ligand-Induced Conformational Changes and Conformational Dynamics in the Solution Structure of the Lactose Repressor Protein. J. Mol. Biol. 2008, 376 (2), 466–481. https://doi.org/10.1016/j.jmb.2007.11.067.

(6) Swint-Kruse, L. Plasticity of Quaternary Structure: Twenty-Two Ways to Form a LacI Dimer. Protein Sci. 2001, 10 (2), 262–276. https://doi.org/10.1110/ps.35801.

(7) Swint-Kruse, L.; Matthews, K. S. Lac Operon. In Encyclopedia of Biological Chemistry; Elsevier, 2013; pp 694–700. https://doi.org/10.1016/B978-0-12-378630-2.00249-8.

(8) Tan, J. S.; Ramanan, R. N.; Ling, T. C.; Mustafa, S.; Ariff, A. B. The Role of Lac Operon and Lac Repressor in the Induction Using Lactose for the Expression of Periplasmic Human Interferon-Α2b by Escherichia Coli. Ann. Microbiol. 2012, 62 (4), 1427–1435. https://doi.org/10.1007/s13213-011-0394-3.

(9) García, Y. Modelación Matemática de La Variabilidad Fenotípica Del Sistema de Regulación Lac En Escherichia Coli, Centro de Investigación en Matemáticas (CIMAT), 2012. https://cimat.repositorioinstitucional.mx/jspui/handle/1008/256.

(10) Munro, P. D.; Ackers, G. K.; Shearwin, K. E. Aspects of Protein–DNA Interactions: A Review of Quantitative Thermodynamic Theory for Modelling Synthetic Circuits Utilising LacI and CI Repressors, IPTG and the Reporter Gene lacZ. Biophys. Rev. 2016, 8 (4), 331–345. https://doi.org/10.1007/s12551-016-0231-9.

(11) Daber, R.; Sharp, K.; Lewis, M. One Is Not Enough. J. Mol. Biol. 2009, 392 (5), 1133–1144. https://doi.org/10.1016/j.jmb.2009.07.050.

(12) Grossman, T. H.; Kawasaki, E. S.; Punreddy, S. R.; Osburne, M. S. Spontaneous cAMP-Dependent Derepression of Gene Expression in Stationary Phase Plays a Role in Recombinant Expression Instability. Gene 1998, 209 (1–2), 95–103. https://doi.org/10.1016/S0378-1119(98)00020-1.

(13) Terpe, K. Overview of Bacterial Expression Systems for Heterologous Protein Production: From Molecular and Biochemical Fundamentals to Commercial Systems. Appl. Microbiol. Biotechnol. 2006, 72 (2), 211–222. https://doi.org/10.1007/s00253-006-0465-8.

(14) Semsey, S.; Jauffred, L.; Csiszovszki, Z.; Erdossy, J.; Steger, V.; Hansen, S.; Krishna, S. The Effect of LacI Autoregulation on the Performance of the Lactose Utilization System in Escherichia Coli. Nucleic Acids Res. 2013, 41 (13), 6381–6390. https://doi.org/10.1093/nar/gkt351.

(15) Andrecut, M.; Cloud, D.; Kauffman, S. A. Monte Carlo Simulation of a Simple Gene Network Yields New Evolutionary Insights. J. Theor. Biol. 2008, 250 (3), 468–474. https://doi.org/10.1016/j.jtbi.2007.10.035.

(16) Wierling, C.; Kühn, A.; Hache, H.; Daskalaki, A.; Maschke-Dutz, E.; Peycheva, S.; Li, J.; Herwig, R.; Lehrach, H. Prediction in the Face of Uncertainty: A Monte Carlo-Based Approach for Systems Biology of Cancer Treatment. Mutat. Res. Toxicol. Environ. Mutagen. 2012, 746 (2), 163–170. https://doi.org/10.1016/j.mrgentox.2012.01.005.

(17) Leigh, J. W.; Bryant, D. Monte Carlo Strategies for Selecting Parameter Values in Simulation Experiments. Syst. Biol. 2015, 64 (5), 741–751. https://doi.org/10.1093/sysbio/syv030.

(18) Siu, N. A Monte Carlo Method for Multiple Parameter Estimation in the Presence of Uncertain Data. Reliab. Eng. Syst. Saf. 1990, 28 (1), 59–98. https://doi.org/10.1016/0951-8320(90)90043-M.

(19) Vinnakota, K. C.; Bugenhagen, S. M. Optimization and Parameter Estimation, Genetic Algorithms. In Encyclopedia of Systems Biology; Dubitzky, W., Wolkenhauer, O., Cho, K.-H., Yokota, H., Eds.; Springer New York: New York, NY, 2013; pp 1600–1604. https://doi.org/10.1007/978-1-4419-9863-7_291.

(20) Schmelcher, M.; Powell, A. M.; Camp, M. J.; Pohl, C. S.; Donovan, D. M. Synergistic Streptococcal Phage λSA2 and B30 Endolysins Kill Streptococci in Cow Milk and in a Mouse Model of Mastitis. Appl. Microbiol. Biotechnol. 2015, 99 (20), 8475–8486. https://doi.org/10.1007/s00253-015-6579-0.

(21) Abdel-Rahman, M. A.; Hassan, S. E. D.; El-Din, M. N.; Azab, M. S.; El-Belely, E. F.; Alrefaey, H. M. A.; Elsakhawy, T. One-Factor-at-a-Time and Response Surface Statistical Designs for Improved Lactic Acid Production from Beet Molasses by Enterococcus Hirae Ds10. SN Appl. Sci. 2020, 2 (4), 573. https://doi.org/10.1007/s42452-020-2351-x.

(22) Geiger, E. O. Statistical Methods for Fermentation Optimization. In Fermentation and Biochemical Engineering Handbook; Elsevier, 2014; pp 415–422. https://doi.org/10.1016/B978-1-4557-2553-3.00021-0.

(23) Ashyraliyev, M.; Fomekong-Nanfack, Y.; Kaandorp, J. A.; Blom, J. G. Systems Biology: Parameter Estimation for Biochemical Models: Parameter Estimation in Systems Biology. FEBS J. 2009, 276 (4), 886–902. https://doi.org/10.1111/j.1742-4658.2008.06844.x.

Results

Model building

Monte Carlo Simulation

Response Surface Model

References

First model

An analysis of the behavior exhibited by each species within the IPTG induction model was conducted, aiming to align their behavior with literature-derived expectations. Figure 1 visually represents how the IPTG concentration diminishes in the medium as it permeates the bacterial cells. Notably, monomers and dimers of the repressor associate, forming tetramers, leading to a predominant concentration of tetramers for a considerable duration. Around the 4-hour mark, we observe an increase on the endolysin production, this occurs as a consequence of the deactivation of the tetramer repressor by forming the tetramer-IPTG complex. At the 5 hour mark the protein production has reached its maximum point. This observed behavior concurs with the laboratory experiments. After being stable for a few moments it begins to decrease the endolysin concentration.
Additionally, this behavior in Figure 1 underscores the critical role of tetramer formation in influencing endolysin production within the model. It substantiates the necessity to consider these dynamic shifts in the repressor states to understand and predict the endolysin concentration dynamics.

<<<<<<< HEAD

=======

Figure 1. Visualization of each species' behavior in the induction model. Accomplished using the estimated parameters for LysCSA13.

>>>>>>> 74d7c07d (Cambios)

As outlined in the methodology, we encountered a limitation regarding certain parameters that were not available in existing literature. This scarcity had a significant impact on the system, prompting us to emphasize the effect of minor variations in kinetic parameters. To achieve this, we employed a Monte Carlo approach and generated graphs reflecting a vast solution space. This approach aimed to delineate the operational ranges within which the model adheres to its assumptions, enabling accurate predictions. In total, 2000 distinct scenarios were generated, each incorporating a standard deviation of 0.05, and subsequently plotted for analysis.

Figure 2.Montecarlo mechanism of LysK induction with 200 nM of IPTG and standard deviation of 0.05.

In our Monte Carlo simulations, we observed that the anticipated concentration of LysK (figure 2) falls within the range of 8000-15000 nM, while for LysCSA13 (figure 3), it lies between 1500-3000 nM. Intriguingly, our simulations unveiled a captivating insight—the potential protein yield of the final endolysin surpasses that of the other species by a staggering factor. By considering the biological variability and comprehending how different factors within the cell influence the protein yield, we can foresee the conditions for maximum protein production. This not only enhances our understanding of the system but also provides valuable insights for scaling up production by predicting the potential protein yield under different conditions. This revelation not only underscores the significance of considering interactions within the system but also opens up exciting possibilities for optimizing protein yield in our bactericidal system.

Figure 3.Montecarlo mechanism of LysCSA13 induction with 200 nM of IPTG and standard deviation of 0.05.

To enhance the confidence and accuracy of the model, our approach involved conducting a focused set of experiments. After confirming the satisfactory performance and alignment of the model with expected behavior, we proceeded with parameter estimation utilizing the deterministic IPTG induction model. Specifically, protein induction experiments at different IPTG concentrations (e.g., 0.2, 0.5, and 1 mM) combined with data obtained from SDS-PAGE gel electrophoresis analyzed using ImageJ served as inputs for the genetic algorithm (GA) to estimate the missing parameters previously mentioned. Table 1 presents the outcomes of this optimization process. The objective was to identify the specific model that optimally fits the data for various concentrations.

Table 1.Estimated kinetic constants for LysK and LysCSA13 at 200 nM IPTG induction. Determined using a GA.

Endolysin	k6	k7	k10	k11	k16	k17	k18	k19
LysK	172806.1332	0.004646925019	8289.871938	1.401669318	7.36E-07	2988.006153	9.87E-08	2411.253642
LysCSA13	528929.5329	0.004682362659	5870.331195	1.596701205	3.94E-07	2284.487239	1.00E-07	1960.781346
Endolysin	k20	k21	k22	k23	k27	k29	k30	k31
LysK	0.1361368059	0.07776261945	0.4241461762	13.7090742	0.08567815958	0.05925751146	0.05214786934	0.8233674927
LysCSA13	0.1283765476	0.2261425676	0.3737824434	10.93457685	0.1271362511	0.1315557746	0.1212103757	1.531937249

Second Model

The results section presents an analysis of a three-level full factorial response surface model investigating the multifaceted effects of varying protein concentrations of LysK, LysCSA13, and LysSS on bactericidal activity towards S. aureus, S. agalactie, and E. coli. The experiment explored all possible combinations of these factors, aiming to elucidate the relationship between protein concentrations and their responses against the targeted bacterial strains.

Figure 5. Response surface graph of coded values, from -1 to 1, representing the low and high limits of protein concentration for LysK and LysCSA13, which are 875-2625 μg/mL and 500-1500 μg/mL, respectively. The response is quantified in terms of the percentage of S. aureus killed.

The quadratic model for bactericidal activity in S. aurers after treatment is significant, with LysK and LysCSA13 being significant factors affecting bactericidal activity. Both LysK and LysCSA13 contribute positively to bactericidal activity, with their quadratic effect being particularly influential for LysCSA13. The impact of LysSS is not visually represented as a distinct axis on the graph. However, there is a control mechanism to adjust this graph at a specific LysSS level. Despite this control, altering the LysSS coded values does not significantly influence the graph, and thus it is not shown. This lack of influence is attributed to the removal of LysSS variables from the model, justified by their observed statistical insignificance in our analysis.

Figure 6. Response surface graph of coded values, from -1 to 1, representing the low and high limits of protein concentration for LysK and LysSS, which are 875-2625 μg/mL and 25-75 μg/mL, respectively. The response is quantified in terms of the percentage of S. agalactie killed.

The 2FI model for bactericidal activity is significant, with LysK (A), LysSS (C), and the interaction term AC being significant factors affecting bactericidal activity. Interestingly, the model establishes that interaction between LysK and LysSS has a synergistic effect on bactericidal activity.

This model aligns with reality, considering that LysK exhibits a broad spectrum of target bacteria, with S. agalactiae being one of them. The potential synergistic effect observed with LysSS could be attributed to its low protein concentration, suggesting that it might contribute to the overall effect observed with LysK.

Figure 7. Response surface graph of coded values, from -1 to 1, representing the low and high limits of protein concentration for LysK and LysCSA13, which are 875-2625 μg/mL and 500-1500 μg/mL, respectively. LysSS is coded at -1 representing 25 μg/mL. The response is quantified in terms of the percentage of E. coli killed.

Figure 8. Response surface graph of coded values, from -1 to 1, representing the low and high limits of protein concentration for LysK and LysCSA13, which are 875-2625 μg/mL and 500-1500 μg/mL, respectively. LysSS is coded at 0.0335 representing ≈ 50 μg/mL. The response is quantified in terms of the percentage of E. coli killed.

Figure 9. Response surface graph of coded values, from -1 to 1, representing the low and high limits of protein concentration for LysK and LysCSA13, which are 875-2625 μg/mL and 500-1500 μg/mL, respectively. LysSS is coded at 1 representing 75 μg/mL. The response is quantified in terms of the percentage of E. coli killed.

Even though the ANOVA analysis indicated the lack of statistical significance for the terms associated with LysSS (C values), the graphical representations above clearly illustrate subtle variations in the response when altering the C value, represented by coded values ranging from -1 to 1 (25 to 75 μg/mL). Intriguingly, the model suggests a bactericidal effect of LysCSA13 on E. coli, contradicting its known target specificity towards S. aureus.

It's important to note that statistical significance and practical significance are not always the same. The ANOVA results may indicate that the term C associated with LysSS is not statistically significant, which means it doesn't contribute significantly to the model's explanatory power based on the available data. However, this doesn't necessarily imply that the C term has no practical or real-world effect on the response. The graphical representation, where changes in C values result in minimal changes in the response, could suggest that the relationship between C and the response is relatively weak but still present.

Optimization

With a comprehensive understanding of the dynamics among the crude extracts of the endolysins for targeting different bacteria, our final goal is to create a highly effective product capable of combating a spectrum of bacterial strains. To achieve an optimal composition, we aim to maximize the bactericidal activity through a general mathematical model equation. This model will yield coefficients that guide the formulation, allowing us to fine-tune and optimize the composition for superior antibacterial effectiveness. In doing so, we aim to develop a potent product that addresses a broad range of bacterial challenges effectively:

Coded values (see table 4 in methods)

Remarkably, the above mathematical model will grant flexibility during product development and process scalability!

Discussion

We observed consistent patterns in the graphs generated during the GA simulation, consistently displaying three lines corresponding to the different induction concentrations. Notably, each graph consistently stabilized at a non-zero point, suggesting that our solution space might be overly complex or potentially suffering from overfitting. It's plausible that we have an excess of parameters or there could be some repetition in the parameters. The persistent stabilization at a non-zero point raised questions about whether we were stuck in a local error¹ or dealing with a highly intricate system, with the latter being the more probable scenario. However, this was not a major concern as the coefficients effectively described protein production under varying conditions.

Our initial parameter configuration, as depicted in Figure 2, revealed suboptimal results. This could be attributed to the intricate nature of the equation system, potentially considering interactions that override each other. The complexity of the equation system implies a highly nonlinear solution space characterized by numerous local minima and maxima.^2–4 Furthermore, due to time constraints our existing set of experimental data proved to be insufficient for comprehensive optimization. Further investigations could encompass diverse pH levels, temperature variations, and a broader range of extracellular IPTG concentrations to enrich the dataset and enhance model accuracy.

Regardless of the distinct graphs we observed, it's noteworthy that the parameter estimation consistently converged to a specific point. Across some of the graphs, several generations of iterations appeared as though the optimization did not show apparent improvement. This observation hinted at the possibility that the coordinates leading to the optimum for each line, specifically the minimum value, remained consistent. The convergence to the same point in parameter estimation suggested a stable and optimal solution for each line, indicating a reliable estimation process.⁵

In the context of the response surface for E. coli, we purposely integrated the C variables (LysSS) into the model during the experimental design. This decision was motivated by the expectation that LysSS, tailored to target and lyse E. coli, would exert a noteworthy influence on the response variable, namely the bactericidal activity (%) against E. coli. However, contrary to our assumption, the statistical analysis revealed that the C variables did not possess significant impact within the model.

One plausible explanation for the lack of statistical significance of LysSS could be attributed to limitations in the experimental setup. In particular, the observed low significance might be linked to the limited availability of total protein concentration during the experimental procedures in the lab. We suspected that the endolysin might possess cytotoxicity, hindering the production of sufficient protein quantities essential for the necessary experiments in our model. Due to the extremely low total protein concentration obtained, we had to work with highly diluted crude extract concentrations, resulting in an exceedingly low endolysin concentration within the sample. With very low concentrations due to limitations mentioned, the impact of the endolysin on E. coli could be negligible, and hence, statistically non-significant.

References

(1) Ashyraliyev, M.; Fomekong-Nanfack, Y.; Kaandorp, J. A.; Blom, J. G. Systems Biology: Parameter Estimation for Biochemical Models: Parameter Estimation in Systems Biology. FEBS J. 2009, 276 (4), 886–902. https://doi.org/10.1111/j.1742-4658.2008.06844.x.

(2) Fernández-Castané, A.; Caminal, G.; López-Santín, J. Direct Measurements of IPTG Enable Analysis of the Induction Behavior of E. Coli in High Cell Density Cultures. Microb. Cell Factories 2012, 11 (1), 58. https://doi.org/10.1186/1475-2859-11-58.

(3) Marbach, A.; Bettenbrock, K. Lac Operon Induction in Escherichia Coli: Systematic Comparison of IPTG and TMG Induction and Influence of the Transacetylase LacA. J. Biotechnol. 2012, 157 (1), 82–88. https://doi.org/10.1016/j.jbiotec.2011.10.009.

(4) Oehler, S. Induction of the Lac Promoter in the Absence of DNA Loops and the Stoichiometry of Induction. Nucleic Acids Res. 2006, 34 (2), 606–612. https://doi.org/10.1093/nar/gkj453.

(5) Marseguerra, M.; Zio, E.; Podofillini, L. Model Parameters Estimation and Sensitivity by Genetic Algorithms. Ann. Nucl. Energy 2003, 30 (14), 1437–1456. https://doi.org/10.1016/S0306-4549(03)00083-5.