Overview

To gain a deeper understanding of the specific mechanisms of various circuits in this project, verify the feasibility of the project, predict outcomes and collaborate with the wet lab team for mutual improvement, we established a series of mathematical models, including a model for the infection of Sinorhizobium in legume roots and the formation of nodules, a nitrogen regulatory protein-DNA docking model, a nitrogen diffusion model, two models for predicting DHA synthesis and a suicide circuit model.

Our modeling work is closely tied to the project. The Sinorhizobium infection model simulates the formation of nodule. The nitrogen diffusion model and DHA synthesis model calculate the final DHA production yield. The protein-DNA docking model pertains to the regulation of the DHA synthesis circuit and confirms the broad applicability of our project. The suicide circuit model ensures the project's biosafety by verifying the effectiveness of the suicide mechanism.

Figure 1 Overview of dry lab

Nodule Formation Model

Background

We aim to have the Sinorhizobium infect the soybean root system to establish an anaerobic synthesis factory. To better calculate the DHA yield from the nodules, we need to estimate the number of Sinorhizobium present within the nodules. Existing studies mainly focus on describing the infection process and mechanism of Sinorhizobium in a qualitative manner. However, to calculate DHA production, we require a method to quantitatively determine the number of Sinorhizobium and attempt to simulate nodule formation. Literature review indicates a standardized process for nodule formation. Given the structured infection process, we opted to build the model using a 3D cellular automaton. By iterating over time, we obtained quantitative data on the growth of nodules and the number of Sinorhizobium at various stages of the process.

Method & Model

To reasonably simulate the formation of nodules, we first examined the process of how Sinorhizobium infects plants’ roots and later forms nodules. Sinorhizobium initially exist as free-living bacteria in the soil. They are attracted to the plant roots via chemotaxis, establishing a bidirectional interaction with the leguminous host plant. The bacteria colonize the root, some attach to root hairs, initiating the plant infection process. Sinorhizobium enter infection threads generated by the plant, moving downward until they encounter developing nodules. The bacteria are then engulfed by plant cells, differentiating into nitrogen-fixing bacteroids. Finally, as the nodule ages, Sinorhizobium are released back into the soil to complete their life cycle.

Figure 2 Formation of the nodules

Given the distinct stages of nodule formation, we recognized the suitability of using a 3D cellular automaton (CA) to represent these states and establish rules for state transitions, completing the model construction.

A cellular automaton (CA) is a grid-based mathematical model used to simulate global behavior emerging from local interactions. Each cell occupies a position on the grid and can be in one of a finite number of states. The state of each cell evolves over time based on a set of simple local rules, determined by the cell's current state and the states of its neighboring cells. Despite the simplicity of these rules, cellular automata can exhibit complex behavior.

A 3D CA extends the concept into three dimensions, allowing for the simulation of systems with 3D structures. In this case, neighbors can be adjacent in all six directions—up, down, left, right, front, and back—making it more precise for simulating phenomena in three-dimensional space.

The formation of nodules is a complex 3D process involving bacterial movement, root growth, infection thread formation, and cell differentiation. A 3D CA is well-suited to simulate this process because the interactions between free-living Sinorhizobium in the soil and the plant root system occur in three dimensions. By stimulating the chemotaxis-driven bacterial movement, infection path formation, and nodule growth in a 3D grid, we can accurately represent these processes. Different cell and environmental states can be simulated using distinct cellular states and hierarchical structures, such as free-living bacteria, bacteria attached to the roots, bacteria that infects root as thread, and N2 fixed bacteria. Each state interacts with others through simple but exacting rules, allowing for an overall representation of the complex biological processes. The dynamic nature of nodule formation is inherently temporal, and the 3D CA models the time-stepped evolution of bacteria, from free-living cells to nitrogen-fixing bacteroids and mature nodules. Therefore, after considering other methods such as kinetics, we selected the cellular automaton to build this model.

Model Assumptions

Spatial Assumptions:
- 3D Grid: A 3D discrete grid represents the space where root and bacterial growth occurs. Each cell in the grid represents a small region of the system and can be in different states (e.g., empty, root, bacteria).
- Boundary Conditions: The grid has finite boundaries, meaning all growth and expansion are confined within the grid.
Neighborhood Rules:
- Moore Neighborhood: Each cell's state evolution depends on its neighbors. The 3D Moore neighborhood assumes that each cell has 26 neighboring cells, including those in front, behind, and in depth (forming a 3×3×3 cube with the center removed).
Time Assumptions:
- Discrete Time Steps: The model progresses in discrete time steps. At each time step, all cells are updated based on their current states and the states of their neighbors.
Root Structure Assumptions:
- Fixed Root Shape: The root is modeled as a slender cylindrical structure growing from the center of the grid. Its shape remains unchanged during the simulation, providing attachment points for the bacteria.
Self-Organizing Behavior:
- Emergence of Global Behavior: Through local rules and random effects, the CA evolves to exhibit global behaviors such as nodule formation. This complexity emerging from local interactions is a key feature of CA models.
Nodule Process Assumptions:
- Ideal Conditions: Under ideal conditions, the growth of nodule strictly follows the mechanisms described in the literature.

Basic States

We strictly followed the literature to define the states, with the only simplification of merging differentiated and undifferentiated nodule cells:

EMPTY (0): An empty cell, representing a location with no activity in the grid.
ROOT (1): Root cells representing the plant’s root system
BACTERIA (2): Free-living Sinorhizobium, representing bacteria not yet attached to the root.
ATTACHED (3): Attached Sinorhizobium, representing those adhering to root hairs.
INFECTION_THREAD (4): Infection threads through which Sinorhizobium enter the root.
NODULE (5): Sinorhizobium in mature nodules, either undifferentiated or differentiated.
N2_FIXING (6): Sinorhizobium in mature nodules that can fix nitrogen, representing fully mature cells capable of nitrogen fixation.

Cellular Automata Rules

Sinorhizobium-Root Interaction:
- Initially, the grid contains root cells (state 1) and free-living Sinorhizobium randomly(state 2).
- Rule 1: If Sinorhizobium (state 2) are adjacent to a root cell (state 1) in the Moore neighborhood, there is a probability of being attracted to the root and changing state to attached Sinorhizobium (state 3).
The Formation and Growth of Infection Thread:
- Rule 2: Once Sinorhizobium attach to roots (state 3), there is a probability of generating an infection thread (state 4), growing inward within the root cells.
- Rule 3: The infection thread (state 4) grows along the root hair until it reaches the nodule primordium.
Nodule Formation and Differentiation:
- Rule 4: Upon reaching a suitable location, there is a probability of the transition from infection thread’s cells to nodule cells (state 5), initiating nodule development.
- Rule 5: Once the nodule cells contact the infection thread, there is a probability of differentiating and enter the nitrogen-fixing stage (state 6).
Nitrogen Fixation and Sinorhizobium Release:
- Rule 6: When the surrounding oxygen concentration decreases, there is a probability of the nitrogen fixation (state 6) begins.
- Rule 7: After a certain number of time steps, there is a probability of the Sinorhizobium releases back into the soil (state 0), completing the life cycle.

Results

We programmed the CA model using Python and obtained the following results.

Figure 3 Nodule formation result1

If the initial amount of Sinorhizobium is limited, the entire root undergoes a process whereby the Sinorhizobium initially infects a small area and then gradually spreads and forms lines of infection that eventually form nodules. Over time, nitrogen-fixing nodules appear, which is consistent with actual nodule growth patterns.

Additionally, we tracked the number of cells in each state as a function of iterations, as shown in the graph below:

Figure 4 Nodule component count over iterations 1

We hypothesize that the nitrogen-fixing bacteria (state 6, represented in brown) are likely responsible for DHA production under nitrogen regulation, making them the focal point of our study.

By examining the growth curves, we observe that uninfected root cells gradually decrease, and free-living bacteria only appear in the initial stages near the root system. After contacting the root, they quickly initiate infection. During the infection burst, the number of infection threads stabilizes as the root decreases and nodules increase. After a period of time, the Sinorhizobium reaches saturation, converting into nitrogen-fixing bacteria that can initiate DHA synthesis. In the later stages of growth, the number of infection threads declines, and the nodule matures, consisting primarily of nodule bacteria and Sinorhizobium, matching the actual development of nodules.

We also explored the effect of the initial Sinorhizobium number on nodule formation by increasing the Sinorhizobium count tenfold：

Figure 5 Nodule formation result 2

The results show that the speed of DHA-producing Sinorhizobium formation significantly increased, indicating that the number of Sinorhizobia introduced is crucial for determining how quickly nitrogen-fixing and DHA-producing Sinorhizobium mature.

The curve below shows how each state of Sinorhizobium changes over iterations:

Figure 6 Nodule component count over iterations 2

We observed that the time required for infection was significantly shortened, and the plateau phase of infection threads disappeared. The rapid growth of nodules broke the balance during the mid-infection phase, leading to a sharp increase in the number of DHA-producing Sinorhizobium.

Contribution

Stimulating the infection process of Sinorhizobium in plants to form nodules can help our team better understand the mechanisms of nodule formation and further quantify DHA synthesis. Current research offers little in the way of quantitative calculations regarding nodules, and we have established a precedent for computational modeling in this area. Future teams interested in Sinorhizobium can refer to our model and make improvements.

Discussion

We initially wanted to use 2D cellular automata, but found that it could not reflect the uniqueness of the three-dimensional structure of the rhizoma, so we switched to 3D. We started with a simplified rule set, through several rounds of iterations, we finally came up with this complete and detailed procedure after referring to a large number of literature. However, our model does have limitations that can be addressed. For instance, the discussion about the formation mechanisms of the spherical structure of nodules is lacking in the literature, which means our model does not achieve visual accuracy in this aspect—which we find regrettable.

Protein-DNA Docking Model

Background

The nitrogen response is the most critical regulatory aspect of the entire project because it initiates the ultimate goal of the nodule factory: to synthesize DHA before soybean harvest. Our team needs to identify the nitrogen-responsive components of the circuit that sense changes of nitrogen content in the soil, which initiates the deletion of phaC2 and the expression of the DHA synthesis enzyme cluster.

Through literature review, we found that in Rhizobium etli, the NtrC (Nitrogen regulatory protein C) plays a significant role as a transcription factor in the nitrogen regulatory system. It typically regulates the expression of downstream genes by binding to specific DNA sequences, particularly at the binding sites located upstream of the glnK/amtB operon. However, there is limited research on nitrogen regulation in our chassis, S. fredii CCBAU45436, and we seek to determine whether similar regulatory mechanisms exist that can be leveraged for nitrogen regulation designs.

According to the literature, the components of the glnK/amtB operon in Rhizobium etli are as follows:

-24/-12 Region: This segment comprises a 16-base pair sequence located upstream of the glnK/amtB transcription start site, with the sequence TGGCACGTCCTTTGAA, situated between -26 and -11. This sequence shows significant homology to the consensus sequence of σN (σ54) dependent promoters (-24 GG/-12 GC).
NtrC Binding Site: In the glnK/amtB promoter sequence, the binding site for NtrC is located between -76 and -91.
Integration Host Factor (IHF) Binding Site: A potential IHF binding site has also been identified between -37 and -72, which has been shown to affect the expression of the glnHPQ operon in E. coli.

These sites are key elements in regulating the transcription of the glnK/amtB operon, where NtrC must bind to these specific sites to initiate gene expression.

Figure 7 Rhizobium etli glnK and amtB genes and partial tesB gene

Method & Model

To explore the relevant mechanisms, we first searched for similar proteins and regulatory sequences. And then we predicted the structure of the NtrC protein and performed molecular docking predictions to score for the protein-DNA interaction by HDOCK. Finally, we visualized the results with PyMOL. After identifying homologous operons and regulatory proteins, our focus shifted to validating the binding of DNA and proteins to determine if a similar effect as observed in Rhizobium etli can be achieved.

Protein-DNA docking involves predicting how proteins interact with DNA at the atomic level. It is a key technique in bioinformatics and molecular biology for studying protein-DNA complexes, which are crucial for processes such as transcription, replication and DNA repair. This docking simulation aids in determining the binding sites, orientation and affinity between the molecules. The specific engineering process is illustrated in the accompanying diagram.

Pipeline

Identifying Operon Gene Sequences
We first located the complete genome of S. fredii CCBAU45436 on NCBI. Using SnapGene, we searched for the testB gene to find the similar glnK/amtB operon in Rhizobium etli. By comparing the sequences upstream of the gene, we identified the corresponding NtrC binding sites and the core region of the glnK/amtB operon.

Figure 8 S. fredii CCBAU45436 glnK and amtB genes and partial tesB gene
Identifying Homologous NtrC Proteins
We found the NtrC protein of S. fredii CCBAU45436 on NCBI, which is homologous to the one mentioned in the literature. Since there is no structural data for NtrC family proteins in the PDB, we used AlphaFold to predict the structures of the NtrC proteins from both types of Rhizobia, yielding satisfactory confidence levels.
HDOCK Docking
The ultimate goal of this model is to predict protein-DNA binding. This prediction differs from traditional protein-protein surface pairing algorithms, as the binding typically relies on smaller contact surfaces, including specific hydrogen bonds and electrostatic interactions, particularly in the recognition sequences formed in the major and minor grooves of DNA. Predicting protein-DNA binding has high specificity because the DNA sequence and the geometric features of the major and minor grooves are relatively fixed, making it easier to identify binding sites. Using HDOCK, which combines template-based and free docking, allows for efficient prediction of protein-DNA binding, especially useful when we have structural information on the predicted protein binding sites for subsequent comparative analysis.

HDOCK's main components include:
- Hybrid Docking Algorithm: Combines template-based docking with de novo docking (free docking) for high precision.
- FFT-based Global Search: Utilizes the Fast Fourier Transform (FFT) algorithm for a global search of potential docking configurations, ensuring a variety of possible interactions are considered.
- Scoring Function: A complex scoring function evaluates docking configurations by considering shape complementarity, electrostatics and hydrogen bonding.
- Flexibility: HDOCK allows for small flexibilities in protein and DNA structures, enabling more accurate predictions of biologically relevant conformations.
Figure 9 The pipeline of HDOCK

To validate the credibility of this model, we first used HDOCK to verify the NtrC protein and glnK/amtB operon from Rhizobium etli, which have been discovered in the literature. Building on that, we input the AlphaFold predicted protein (.pdb) and operon sequences into HDOCK, running the program on the server to obtain the ten highest probability docking results.
Visualization with PyMOL
By importing the docking results into PyMOL for visualization, we analyzed the binding sites and residue numbers, comparing them to the results from SnapGene. Besides, we also topologized the protein and DNA to prepare for further molecular dynamics simulations.

Results

AlphaFold Prediction Results:
Prediction for Rhizobium etli:

Figure 10 Predicted structure for Rhizobium etli

Prediction for S. fredii CCBAU45436:

Figure 11 Predicted for S. fredii CCBAU45436

The protein sequences show that the NtrC proteins from Rhizobium etli and S. fredii CCBAU45436 consist of 484 and 483 amino acid residues, respectively. The structural predictions yielded confidence levels of 0.73 and 0.74, indicating good reliability. Observations confirm they are homologous proteins with very similar structures. Additionally, visualization suggests the protein has a palm-like shape with a longitudinal groove, hypothesized to be the binding site for the protein as a transcription factor with DNA, though further validation is required.
HDOCK Prediction Results:
Rhizobium etli:

Figure 12 The HDOCK result of Rhizobium etli

S. fredii CCBAU45436:

Figure 13 The HDOCK result of S. fredii CCBAU45436

HDOCK generates multiple scenarios; here we list the top ten best-performing results. The docking scores calculated by HDOCK are based on knowledge-driven iterative scoring functions, ITScorePP or ITScorePR. The docking score reflects the likelihood of the predicted binding model: more negative values indicate more stable binding between the two molecules. For protein-protein or protein-nucleic acid complexes, scores of -200 or better represent high confidence. In the previously validated model for Rhizobium etli, the docking score was -201.42, suggesting that DNA and protein can bind well, which was consistent with the literature and validated the model for further predictions. For S. fredii CCBAU45436, the docking score was -222.99, indicating even better binding.

The confidence score, based on the docking score, quantifies the likelihood of binding between the two molecules. It is calculated using the formula:

'$$Confidencescore = 1.0/[1.0 + e^{(0.02 * (Docking_Score + 150))}]$$

A confidence score exceeding 0.7 indicates a high likelihood of binding; scores between 0.5 and 0.7 indicate potential binding; while scores below 0.5 suggest low binding likelihood. In Rhizobium etli, the confidence score for binding was 0.7366 (>0.7), indicating that the model is consistent with the literature and can continue to predict. In S. fredii CCABU45436, the confidence score for binding was 0.8115, indicating very good binding results.

Both metrics suggest that the model has good predictive capability, and NtrC can bind to the glnK/amtB operon in both Rhizobium etli and S. fredii CCBAU45436.
PyMOL Visualization Results:
We imported the predicted binding results generated by HDOCK into PyMOL. Combining this with the residues provided by HDOCK, we displayed the amino acid residues, base indices and names, along with the connecting hydrogen bonds and their lengths at the docking interface.

Rhizobium etli:

Figure 14 Docking of NtrC protein and DNA in Rhizobium etli

S. fredii CCBAU45436:

Figure 15 Docking of NtrC protein and DNA in S. fredii CCBAU45436

From the figures, the NtrC proteins in both bacteria can bind to the glnK/amtB operon, with binding occurring at the grooves as predicted. Additionally, we found that the NtrC protein in Rhizobium etli has two concentrated binding sites with DNA, consistent with the NtrC binding site and the IEF binding site mentioned in the literature. The distal contact site may correspond to the core functional region of the operon. Based on the verification of the literature, our chassis binding model is particularly significant. The base numbers of the DNA bases at the binding interface, as shown in the figure, are almost identical to those predicted in the SnapGene, aligning with the NtrC binding site. However, the distal binding is notably more extensive than that of Rhizobium etli, suggesting a tighter interaction, which may also imply other functions. This aligns with the higher scores from HDOCK, indicating that NtrC in S. fredii CCBAU45436 is not only structurally homologous but also functionally consistent with the described homology in the literature.

Contribution

This model validates the feasibility of the circuit supported by existing literature, which shows that the NtrC protein in Rhizobium etli can bind to and regulate the glnK/amtB operon. Subsequently, we used this circuit to predict whether S. fredii CCBAU45436 possesses a similar mechanism. Ultimately, we identified proteins and DNA with high structural and functional homology, confirming the existence of a nitrogen regulatory transcription factor. This model assists the team in establishing NtrC protein as the most critical nitrogen response regulatory element, significantly accelerating the progress of wet lab and enhancing their confidence.

Discussion

We initially chose many molecular docking models to predict binding, but found that most of the existing models focus on protein-protein docking. After reviewing a large amount of literature and comparing multiple models, we chose the HDOCK model. Furthermore, this model illustrates the broad applicability of our project, indicating that our nodule factory is not limited to the selected chassis bacterium. Through future analysis of other rhizobia, we can reveal that similar regulatory mechanisms are present across Sinorhizobia, thereby continuously expanding our nodule factory.

DHA Yield Prediction Model

Overview

The DHA yield prediction model is designed to simulate and optimize the production of docosahexaenoic acid (DHA) using Sinorhizobium. This model integrates various biochemical pathways and environmental factors that influence DHA synthesis, enabling accurate predictions of yield under different conditions.

When building the Model, we first considered the static situation without considering the time change. In order to ensure the accuracy of the data, FBA (Flux Balance Analysis) Model was established. However, it was found that the time change needs to be considered, so the ODE model which could correspond with timeline was established. After the sensitivity analysis of the ODE model, it was found that the activity of the PUFA synthase has a significant effect on the results and is regulated by Nitrogen, so a Nitrogen Absorption Model is established.

Background

We developed a novel biosynthesis model that utilizes Sinorhizobium for the efficient production of docosahexaenoic acid (DHA). The core of this model is to introduce the pfa biosynthetic gene cluster for increasing the yield of EPA and converting eicosapentaenoic acid (EPA) to DHA, allowing it to produce DHA. Additionally, we deleted the acetyl-CoA pathway that generates polyhydroxybutyrate (PHB), maximizing carbon flow.

In practical experimental operations, validating DHA yield faces challenges due to the long plant growth duration and the difficulty of conducting experiments. Therefore, we predict DHA yield through two kinds of models. Furthermore, since the production rate of DHA is regulated by nitrogen concentration, our model also simulates nitrogen diffusion and the absorption of nitrogen by plant roots to model the process of nitrogen regulation on DHA production.

Method & Model

Nitrogen Absorption Model

To accurately describe the impact of nitrogen on DHA yield, a nitrogen absorption model of Sinorhizobium has been incorporated, as the synthesis rate of PUFA synthase is positively correlated with nitrogen concentration.

Method

Numerical realization of convection and diffusion in soil:
- Diffusion part: the diffusion of nitrogen in soil is calculated using the Laplace operator. The nitrogen concentration matrix is shifted in three directions (x, y, z) to calculate the concentration change of six adjacent points, simulating three-dimensional diffusion.
- Convection part: the convection term is only calculated along the z direction, simulating the migration of nitrogen driven by water flow. The nitrogen concentration in the entire soil grid is updated after each time step dt.
Realization of nitrogen absorption by nodules:
- Calculation of nodule surface: the surface of nodules is extracted by performing a convolution on the six neighboring positions in the nodule region.
- Absorption update: the nitrogen concentration on the nodule surface is updated based on the given absorption rate U, simulating the process of nitrogen being absorbed.
Time curve of nitrogen absorption by nodules:
- A red frame is used to enclose the 3x3 region where the nodules are located. In each time step, the total amount of nitrogen absorbed on the nodule surface is calculated and stored. As time progresses, these absorption values are recorded and a graph is drawn showing the relationship between the total amount of nitrogen absorbed by nodules and time.

Equations

Diffusion Model in Porous Media
- Assumption: Soil can be viewed as a porous medium, where nitrogen is transported to plant roots through diffusion and convection.
- Equation: $$ \frac{\partial C}{\partial t} = D \nabla^2 C - \nabla \cdot (v C) + R(C, t) $$ Where:
- $C$ is the nitrogen concentration.
- $D$ is the diffusion coefficient of nitrogen in the soil.
- $v$ is the soil water flow velocity (convection term).
- $R(C, t)$ is the reaction term for nitrogen, including absorption and reaction.
Root Absorption Model
- Assumption: The Sinorhizobium absorbs nitrogen from the soil through the plants' roots.
- Process: $$ U = \frac{U_{max} C}{K_m + C} $$ Where:
- $U$ is the absorption rate.
- $U_{max}$ is the maximum absorption rate.
- $K_m$ is the half-saturation constant.
Root Nodule Diffusion Model
- Assumption: Nitrogen diffuses within the root nodule and is transferred to the Sinorhizobium.
- Equation: $$ J = -D_{root} \nabla C $$ Where:
- $J$ is the diffusion flux.
- $D_{root}$ is the diffusion coefficient within the root nodule.

Flux Balance Analysis

Introduction

Flux Balance Analysis (FBA) is a mathematical approach used to analyze metabolic networks by optimizing the flow of metabolites through a network, based on stoichiometric relationships. It involves constructing a metabolic model from the organism's genome, where each reaction is described by a stoichiometric matrix. FBA does not require detailed kinetic parameters, as it relies on the assumption of steady-state conditions, where the sum of all fluxes entering and leaving a metabolite is zero. By applying linear optimization, FBA predicts the growth rate or other objectives without needing specific reaction rates or enzyme data.

To ensure the accuracy of the model, we utilized the Cobra library in Python for flux balance analysis. This analysis process includes:

Model Construction: Integrate the metabolic network, defining relevant reactions along with their substrates and products.
Constraint Setting: Establish substrate availability and reaction rates based on experimental conditions.
Objective Function: Set the DHA yield as the optimization target, maximizing its synthesis efficiency.
Optimization Solving: Apply linear programming methods to obtain the maximum theoretical yield of DHA.
Result Analysis: Analyze the contribution of each reaction, identify production bottlenecks, and provide a basis for further optimization.

Through this series of steps, we gain insights into the DHA production capacity of Sinorhizobium under specific conditions, offering theoretical support for future experimental designs.

Model Assumptions

Completeness of Metabolic Pathways: The constructed metabolic network model includes all reactions related to DHA synthesis, and the kinetics of these reactions are reliable.
Reaction Rates: Assume that all metabolic reaction rates can be described by mass balance and remain stable during growth, suitable for linear programming assumptions.
Coupling of Cell Growth and Metabolism: The metabolic flux between cell growth and DHA synthesis is balanced, meaning that cell growth does not limit DHA production, and vice versa.
Constant Environmental Conditions: The environmental conditions considered in the model (such as temperature, pH, and oxygen concentration) remain unchanged during the analysis.
Reversibility of Reactions: Assume that the involved reactions can be reversed under certain conditions during biological conversion, and their equilibrium constants are known.
Enzyme Activity: The enzyme activity required for the reactions in the model is sufficient to drive the reactions toward product formation without rate limitations.
Gene Expression Levels: Assume that the expression levels of introduced genes (such as those related to EPA synthesis) in Sinorhizobium are sufficiently high to support efficient DHA production.
Metabolic Flux Balance: Assume that all metabolic fluxes in the model are balanced, meaning inputs and outputs are equal under steady-state conditions.

Model Establishment

Identifying and Removing PHB Production Reactions:
We used a metabolism model of S. fredii CCBAU45436 according to the metabolic measurements by other researchers (Osbaldo.R, 2007). The model used the following code to delete reactions related to PHB production.
```
phb_reactions = [rxn for rxn in main_model.reactions if 'PHB' in rxn.id or 'PHB' in rxn.reaction]
    for rxn in phb_reactions:
        print(f"Removing Reaction ID: {rxn.id}")
        print(f"Reaction Equation: {rxn.reaction}")
        print()
    main_model.remove_reactions(phb_reactions)
    print("PHB production pathways have been removed from the main model.")
    
```
The following two reactions were removed:
- Removing Reaction ID: PHBS
  Reaction Equation: $ \text{3HBCoA}_{\text{R,c}} \rightarrow \text{CoA}_{\text{c}} + \text{PHB}_{\text{c}} $
- Removing Reaction ID:$PHB_{deg_{1}}$
  Reaction Equation: $ \text{H}_2\text{O}_{\text{c}} + 2.0 \text{PHB}_{\text{c}} \rightarrow \text{3HBB}_{\text{RR,c}} + \text{H}_{\text{c}} $

We referred to the complete metabolic model of Escherichia DH5α (Jonathan M Monk, 2016) to add reactions related to DHA. The following were retrieved:

Reaction ID: DHAD1
Reaction Equation: $ \text{23DHBM}_{\text{c}} \rightarrow \text{3MOB}_{\text{c}} + \text{H}_2\text{O}_{\text{c}} $
Reaction ID: DHAD2
Reaction Equation: $ \text{23DHMP}_{\text{c}} \rightarrow \text{3MOP}_{\text{c}} + \text{H}_2\text{O}_{\text{c}} $
Reaction ID: DHAPT
Reaction Equation: $ \text{DHA}_{\text{c}} + \text{PEP}_{\text{c}} \rightarrow \text{DHAP}_{\text{c}} + \text{PYR}_{\text{c}} $
Reaction ID: DHAt
Reaction Equation: $ \text{DHA}_{\text{e}} \leftrightarrow \text{DHA}_{\text{c}} $
Reaction ID: GMHEPAT
Reaction Equation: $ \text{ATP}_{\text{c}} + \text{GMHEP1P}_{\text{c}} + \text{H}_{\text{c}} \rightarrow \text{ADPHEP}_{\text{DASH,DD,c}} + \text{PPi}_{\text{c}} $

Based on this, we constructed the metabolic flux model for flux balance analysis, and we used Cytoscape for visualization, generating related metabolic network diagrams.

Figure 16 metabolic network diagrams of the model

Differential Equation Model

This model describes the changes in the concentrations of metabolic products over time using differential equations. The synthesis of DHA is proportional to the concentration of EPA, while the secretion rate of DHA is proportional to its intracellular concentration, capped by a maximum value. This maximum value is determined by the maximum flow rate obtained from the Flux Balance Analysis Model, ensuring that we account for temporal changes based on precise values.

Model Assumptions

Dynamic Reactions: The reaction rates are assumed to vary over time, reflecting real-time changes in concentration.
Linear Relationship: It is assumed that the relationship between reaction rates and substrate concentrations is linear, suitable for Michaelis-Menten kinetics.
Initial State:The system is assumed to be in a stable initial state at the start of the reaction, with all component concentrations at equilibrium.
Stable Enzyme Activity: It is assumed that the activity of the enzymes involved in the reactions remains stable throughout the reaction, influenced only by the rates controlled within the model.

Parameter Settings

Parameters	Description	Value	Unit
$Ac_{in}$	Input Concentration of Acetyl-CoA	$1$	$\text{mol/L}$
$Ac_{out}$	Output Concentration of Acetyl-CoA	$1$	$\text{mol/L}$
$diff$	Diffusion Rate	$0.1$	$\text{mol/(L·h)}$
$v_{PHBmax}$	Maximum PHB Production Rate	$1$	$\text{mol/(L·h)}$
$K_{mPHB}$	PHB Michaelis Constant	$0.5$	$\text{mol/L}$
$v_{EPAmax}$	Maximum EPA Production Rate	$1$	$\text{mol/(L·h)}$
$K_{mEPA}$	EPA Michaelis Constant	$0.5$	$\text{mol/L}$
$c_{Ac0}$	Initial Acetyl-CoA Concentration	$0.1$	$\text{mol/L}$
$c_{PHB0}$	Initial PHB Concentration	$0.1$	$\text{mol/L}$
$c_{EPA0}$	Initial EPA Concentration	$0.1$	$\text{mol/L}$
$DHA_{c0}$	Initial DHA Concentration	$0.0$	$\text{mol/L}$
$DHA_{e0}$	Initial DHA Secretion Concentration	$0.0$	$\text{mol/L}$
$k$	DHA Synthesis Rate Constant	$0.5$	$\text{mol/(L·h)}$
$K_{DHA}$	DHA Synthesis Rate Threshold	$5.0$	$\text{mol/L}$
$S_{max}$	Maximum DHA Secretion Rate	$0.5$	$\text{mol/(L·h)}$
$K_{S}$	DHA Secretion Michaelis Constant	$2.0$	$\text{mol/L}$

Table 1 Parameter settings of Differential Equation Model

Notably, the model considers the impact of bacterial growth, where the concentration decreases as the volume increases, reflected in the diffusion term.

Differential Equations

In the process of listing the differential equation, we referred to the metabolic diagram below, and simplified the redundant steps of the shaded region to make the model more concise.

Figure 17 Metabolic diagram of the Differential Equation Model

Change in Acetyl-CoA Concentration: $$ \frac{dcAc}{dt} = Ac_{in} - cAc \cdot (Ac_{out} + v_1 + v_2) - diff $$
Change in PHB Concentration: $$ \frac{dcPHB}{dt} = v_1 - diff $$
Change in EPA Concentration: $$ \frac{dcEPA}{dt} = v_2 - diff $$
Change in Intracellular DHA Concentration: $$ \frac{dDHA_c}{dt} = k \cdot cEPA - \frac{DHA_c}{K_{DHA}} $$
Change in Secreted DHA Concentration: $$ \frac{dDHA_e}{dt} = \frac{S_{max} \cdot DHA_c}{K_{S} + DHA_c} $$

Reaction Rate Equations

PHB Synthesis Rate: $$ v_1 = \frac{v_{PHBmax} \cdot cPHB}{K_{mPHB} + cPHB} $$
EPA Synthesis Rate: $$ v_2 = \frac{v_{EPAmax} \cdot cEPA}{K_{mEPA} + cEPA} $$
DHA Synthesis Rate:
$k$: DHA synthesis rate constant
$K_{DHA}$: Threshold for decreasing DHA synthesis rate
$S_{max}$: Maximum DHA secretion rate
$K_{S}$: Michaelis constant for secretion rate concerning DHA concentration

Results

Nitrogen Contact Model Results

Nitrogen Concentration over Time in the Soil

The visualized 2D slice of nitrogen concentration illustrates the distribution of nitrogen in the soil grid at various time steps. The red rectangle indicates the root nodule's location, while the color map represents the nitrogen concentration. Early time steps show a nitrogen concentration peak at the initial fertilizer point outside the nodule, with diffusion gradually spreading nitrogen outward. Over time, nitrogen reaches the nodule's surface, resulting in absorption, which eventually depletes the nitrogen concentration around the nodules.Nitrogen gradually approaches the nodule, indicated by the spread towards the red rectangle (the nodule zone). This progression reflects the transport dynamics governed by diffusion (with a low diffusion coefficient) and convection (slow velocity), leading to a steady accumulation of nitrogen around the nodule.

Figure 18 Time-series plot of nitrogen diffusion

The figure presents a time-series plot of the total nitrogen absorbed by the nodule as a function of time. The amount of nitrogen absorbed follows a bell curve pattern: it increases rapidly as nitrogen reaches the nodule, peaks, and then slowly declines as the available nitrogen decreases in the surrounding soil.

Figure 19 Time-series plot of the total nitrogen

Analysis:

The rapid initial increase corresponds to nitrogen diffusing to the nodule's surface. As nitrogen reaches the nodule, it is absorbed at a steady rate, causing the concentration on the surface to rise.
The peak absorption rate occurs when the maximum nitrogen flux into the nodule is reached. Beyond this point, nitrogen in the surrounding soil becomes depleted due to absorption and the diminishing nitrogen supply.
The slow decline in absorption is explained by the gradual depletion of nitrogen around the nodule, leaving less available for continued absorption.

The results from both plots indicate that nitrogen diffusion from a localized source outside the nodule is initially slow but accelerates as nitrogen approaches the nodule surface. The nodule begins absorbing nitrogen at a growing rate, which peaks and then decreases as the local nitrogen concentration is consumed.

This model illustrates a classic nutrient absorption process: diffusion-convection governs nitrogen movement through the soil, while the nodule acts as a sink, progressively reducing nitrogen availability. Over time, the system reaches a quasi-equilibrium where nitrogen supply in the soil diminishes as more nitrogen is absorbed by the nodule.

Flux Balance Analysis Results

Reaction	Value
Target Value	10.00
PFK	10.00
PGI	10.00
PGK	-20.00
PGM	-20.00
ATPM	20.00
PYK	10.00
D_LACt2	-20.00
ENO	20.00
TPI	10.00
EX_glcD_e	-10.00
EX_h	20.00
EX_lacD_e	20.00
FBA	10.00
GAPD	20.00
GLCpts	10.00
LDH_D	-20.00

Table 2 The result of Flux Balance Analysis

Analysis of Reaction Impact on DHA Production

In our optimization, we set the reaction with ID FBA (DHA secretion pathway) as the maximization objective. The results excluded reactions with negligible flux, revealing that the reactions listed above significantly influence the total metabolic flux when maximizing DHA yield, which is 10.00 mmol/gDW/hr. However, this number is obviously larger than the actual situation, because due to the influence of time, nitrogen, and so on, the actual yield can not reach the maximum.

But this insight suggests that these reactions listed above are critical points for further analysis in practical production processes, providing a foundation for targeted improvements and optimization strategies.

Differential Equation Model Results

Numerical solutions obtained using MATLAB's ode45 method reveal concentration data at different time points. Simulation results show that the substrate EPA for DHA increases at an appropriate rate and stabilizes over time.

The simulation results indicate that the substrate EPA for DHA production increases at an appropriate rate. This growth aligns with our model's predictions, confirming the efficacy of the metabolic pathways involved in DHA synthesis.

Figure 20 Concentration vs. Time in Differential Equation Model

Sensitivity Analysis

A sensitivity analysis of the above parameters was conducted to explore how parameter values affect DHA yield. The summary of the sensitivity analysis results is as follows:

Figure 21 Sensitivity Analysis of Differential Equation Model

$S_{max}$ and $K_S$ have little effect on DHA concentration:This indicates that, within the current range of model parameters, the maximum secretion rate and Michaelis constant for DHA contribute minimally to the final DHA concentration. These parameters could be fixed in practical applications to simplify the model.
$k$ and $K_{DHA}$ significantly affect DHA concentration:This indicates that, within the current range of model parameters, the maximum secretion rate and Michaelis constant for DHA contribute minimally to the final DHA concentration. These parameters could be fixed in practical applications to simplify the model.

Integrated Model Results Analysis

Assuming that DHA synthase activity is proportional to nitrogen concentration, we can qualitatively analyze that nitrogen concentration will first increase rapidly and then decrease slowly, and since k value is a determining factor that has a significant influence on the yield, DHA secretion rate should first increase rapidly and then stabilize.

Contribution & Discussion

Contribution:This model provides a quantitative tool for understanding the metabolic pathways of Sinorhizobium and can help optimize the production of DHA and its related metabolites. Further validation of the model combined with experimental data can provide strong support for industrial applications.
Reflection of Nitrogen Diffusion and Absorption Model:In our nitrogen diffusion and absorption model, using a uniform diffusion coefficient and convection velocity doesn't account for the spatial and temporal variability in real soil environments. Additionally, the fixed grid size for the nodule also doesn't reflect actual growth dynamics. To enhance the model, we could incorporate spatial variability in soil properties, make absorption rates concentration- or time-dependent, and allow for dynamic changes in nodule size. Validation with experimental data would further strengthen the model's predictive capabilities and improve its representation of nitrogen dynamics.
Emphasizing the Importance of Enzyme Activity Measurement:Given the high sensitivity of the k value to the final result, it is recommended to prioritize accurate measurement of DHA synthase activity in experiments. This will aid in better optimizing the model and enhancing DHA yield.

Suicide Circuits Model

Background

Sinorhizobium gathers at the roots of plants, providing good growth conditions for nitrogen fixation. With synthetic biology, we can make Sinorhizobium capable of producing DHA while fixing nitrogen, which provides better economic benefits for people. However, we do not want our modified Sinorhizobium to escape and spread out of the soil, causing adverse effects such as genetic contamination. With these concerns, we designed the suicide circuit. To verify that our suicide circuit is safe and feasible, we constructed a kinetic model of the net toxin expression.

Method & Model

Suicide Mechanism

Figure22 The dynamics of the kill switch construct

We formulated differential equations based on the properties and interactions of genes.

In the circuit, VapC is a toxin, VapB is an antitoxin. VapB can bind to VapC, blocking the toxin of VapC from functioning. Cas12k is a protein without cutting function and can bind to mCherry under the guidance of sgRNA-mCherry, obstructing the expression of VapC ultimately.

nifH is a promoter in response to hypoxia signals, and glnK is a promoter in response to low nitrogen signals, making VapB highly expressed under low oxygen environment as well as sgRNA-mCherry highly expressed under low nitrogen environment.

As a result, the following four conditions can occur:

1.Low oxygen and high nitrogen -- antitoxin of VapB is expressed, while sgRNA-mCherry is not transcribed. Cas12K cannot bind to mCherry without the guidance of sgRNA-mCherry. Therefore, toxin of VapC is expressed, but Sinorhizobium can survive.

2.High oxygen and high nitrogen -- antitoxin of VapB is not expressed and sgRNAmCherry is not transcribe. At this time, Cas12K cannot bind to mCherry due to the absence of guidance of sgRNA-mCherry. Toxin of VapC is expressed, resulting in death of Sinorhizobium.

3.Low oxygen and low nitrogen -- antitoxin of VapB is expressed and sgRNA-mCherry is transcribed. Cas12K binds to mCherry under the guidance of sgRNA-mCherry. Toxin of VapC is not expressed and Sinorhizobium can survive.

4.High oxygen and low nitrogen – neither VapB nor sgRNA-mCherry is expressed or transcribed. At this time, Cas12K binds to mCherry under the guidance of sgRNA-mCherry. Toxin of VapC is not expressed and Sinorhizobium can survive.

In summary, only when nitrogen fertilizer is applied, which means nitrogen content increases, the Sinorhizobium will commit suicide and die, corresponding to the conditions when it escapes from the soil roots thus causing oxygen content increases.

We considered six different chemical substances: sgRNA-mCherry(sg), mRNA VapC(mC), mRNA VapB(mB), mRNA Cas12K(mK), VapC(C), VapB(B), Cas12K(K), the combination of VapC and VapB(CB), whose concentrations are determined by ordinary differential equations.

Model Assumptions

1.Quasi-steady state approximation: Due to a difference in the time scales involved between binding/unbinding of nifH and oxygen, glnK and nitrogen (happening in the order of milliseconds) and transcription/translation (happening in the order of minutes). We assume that the concentration of nitrogen-activated glnK promoter , oxygen-activated nifH promoter remains constant for time, that is, ${d[glnK-N]}{dt} = 0, \frac{d[nifH-O]}{dt} = 0$.

2.The transcription rates: both the toxin and antitoxin mRNAs were assumed to be approximately similar.

3.The initial concentration: all chemical substances is assumed to be zero.The following table shows the descriptions and values of each parameter.

Parameters setting

The following table shows the descriptions and values of each parameter.

Parameters	Description	Value	Unit
$v_{M_C}$	$M_C$ transcription rate	$0.1333$	$s^{-1}$
$v_{M_B}$	$M_B$ transcription rate	$0.1333$	$s^{-1}$
$v_{sg}$	$sg$ transcription rate	$0.1333$	$s^{-1}$
$v_{M_K}$	$M_K$ transcription rate	$0.1333$	$s^{-1}$
$\beta_{C}$	$M_C$ rate of translation	$0.033$	$s^{-1}$
$\beta_{B}$	$M_B$ rate of translation	$0.139$	$s^{-1}$
$\beta_{K}$	$M_K$ rate of translation	$0.139$	$s^{-1}$
$\delta_{M_C}$	$M_C$ rate of decay	$0.00203$	$s^{-1}$
$\delta_{M_B}$	$M_B$ rate of decay	$0.00203$	$s^{-1}$
$\delta_{M_K}$	$M_K$ rate of decay	$0.00203$	$s^{-1}$
$\delta_{sg}$	$sg$ rate of decay	$0.00203$	$s^{-1}$
$d_C$	$C$ rate of decay	$2.8881 \times 10^{-4}$	$s^{-1}$
$d_B$	$B$ rate of decay	$0.00115524$	$s^{-1}$
$d_K$	$K$ rate of decay	$0.00115524$	$s^{-1}$
$F_{KC}$	Cas12k rate of inhibition	$0.2$	-
$a_{CB}$	CB rate of complex formation	$2 \times 10^{-3}$	$M^{-1} s^{-1}$
$\theta_{sgC}$	$sg$ binding C-gene rate	$7.14972 \times 10^{-3}$	$M^{-1} s^{-1}$
$K_D$	$glnK-N$ response index $nifH-O_2$ response index	$800 \times 10^{-6}$	$s^{-1}$

Table 3 Parameters of suicide circuit model

Equations

Low oxygen and low nitrogen:

\[ \frac{d[M_B]}{dt} = \frac{v_{M_B}}{1 + K_D [O_2]} - \delta_{M_B} [M_B] - \beta_B [M_B] \] \[ \frac{d[sg]}{dt} = \frac{v_{sg}}{1 + K_D [N]} - \delta_{sg} [sg] - \theta_{sgC} [sg] \] \[ \frac{d[M_C]}{dt} = v_{M_C} - \delta_{M_C} [M_C] - \theta_{sgC} [sg] - \beta_C [M_C] \] \[ \frac{d[M_K]}{dt} = v_{M_K} - \delta_{M_K} [M_K] - \beta_K [M_K] \] \[ \frac{d[C]}{dt} = \beta_C [M_C] - d_C [C] - \alpha_{CB} [C][B] \] \[ \frac{d[B]}{dt} = \beta_B [M_B] - d_B [B] - a_{CB} [C][B] \] \[ \frac{d[K]}{dt} = \beta_K [M_K] - d_K [K] - F_K C \theta_{sgC} [sg][K] \] \[ \frac{d[CB]}{dt} = a_{CB} [C][B] \]

High oxygen and low nitrogen:

\[ \frac{d[sg]}{dt} = \frac{v_{sg}}{1 + K_D [N]} - \delta_{sg} [sg] - \theta_{sgC} [sg] \] \[ \frac{d[M_C]}{dt} = v_{M_C} - \delta_{M_C} [M_C] - \theta_{sgC} [sg] - \beta_C [M_C] \] \[ \frac{d[M_K]}{dt} = v_{M_K} - \delta_{M_K} [M_K] - \beta_K [M_K] \] \[ \frac{d[C]}{dt} = \beta_C [M_C] - d_C [C] \] \[ \frac{d[K]}{dt} = \beta_K [M_K] - d_K [K] - F_K C \theta_{sgC} [sg][K] \]

Low oxygen high nitrogen:

\[ \frac{d[M_B]}{dt} = \frac{v_{M_B}}{1 + K_D [O_2]} - \delta_{M_B} [M_B] - \beta_B [M_B] \] \[ \frac{d[M_C]}{dt} = v_{M_C} - \delta_{M_C} [M_C] - \beta_C [M_C] \] \[ \frac{d[M_K]}{dt} = v_{M_K} - \delta_{M_K} [M_K] - \beta_K [M_K] \] \[ \frac{d[C]}{dt} = \beta_C [M_C] - d_C [C] - \alpha_{CB} [C][B] \] \[ \frac{d[B]}{dt} = \beta_B [M_B] - d_B [B] - a_{CB} [C][B] \] \[ \frac{d[K]}{dt} = \beta_K [M_K] - d_K [K] \] \[ \frac{d[CB]}{dt} = a_{CB} [C][B] \]

High oxygen high nitrogen:

\[ \frac{d[M_C]}{dt} = v_{M_C} - \delta_{M_C} [M_C] - \beta_C [M_C] \] \[ \frac{d[M_K]}{dt} = v_{M_K} - \delta_{M_K} [M_K] - \beta_K [M_K] \] \[ \frac{d[C]}{dt} = \beta_C [M_C] - d_C [C] \] \[ \frac{d[K]}{dt} = \beta_K [M_K] - d_K [K] \]

Results

We use MATLAB to solve the differential equations, describing the expression of toxins in various environments. The results are as follows:

Figure 23 Result of VapC toxin content under different conditions

It can be seen that only in high nitrogen and high oxygen environment, the net expression of toxins will increase significantly, leading to suicide of Sinorhizobium, in line with expectations.

On this basis, we simulated the whole process of plant growth to further verify the safety and stability of the suicide circuit.

a. At the beginning, Sinorhizobium is introduced: low nitrogen and high oxygen.
b. Development of nodule: Low nitrogen and low oxygen.
c. Application of nitrogen fertilizer: high nitrogen and low oxygen.
d. Sinorhizobium escapes: high nitrogen and low oxygen.

Figure24 Result of VapC toxin in different time periods

The whole process of the toxin content can be observed when the nodule grows in the soil root system. Once the nodule escapes, the toxin content rises quickly, leading to the death of the Sinorhizobium, preventing gene pollution and other problems. It can be suggested that our suicide circuit is safe and feasible.

Sensitivity Analysis

By disturbing the reaction rate constant with different amplitudes, it can be observed that the substances VapC and Cas12k do not respond much to the change of reaction rate, while the production of B and CB has a tendency to decrease and increase by about 3% respectively under the perturbation. It can be seen from the analysis that the disturbance of reaction rate has little effect on the overall effect and does not greatly interfere with the generation of the suicide circuit.

Figure25 Result of sensitivity analysis of parameters

Contribution

This model provides a strong proof for the safety and feasibility of our project, and provides insights for predicting the death of Sinorhizobium, saving the time and energy for experimental verification and providing reference and convenience for subsequent preparation. In addition, we provide a dual-response component model which can be used as an element.

Discussion

At the very beginning of modeling the suicide circuit, we referred to Part: BBa_K3512042. Based on this, we improved their differential equations to make them more close to reality, such as adding a condition that the tendency of newly generated proteins slow down the production of proteins. What’s more, our suicide circuit also involves two response-regulated promoters, the interaction between Cas12k and sgRNA-mCherry, and the mutual binding of VapC and VapB. After the constant modification and improvement of the equations, we finally formulated more complex and rigorous, more parameterized differential equations, which simulates the whole process of the suicide circuit to the maximum extent.

For specific iteration, please read our engineering.

Although this model provides an obvious qualitative prediction of toxin growth in Sinorhizobium, it was unable to achieve further precise quantitative predictions.The selection of coefficient values is largely based on literature estimation, which is inaccurate.

In addition, in low nitrogen and high oxygen environment, it can be seen that although the toxin content tends to stabilize quickly, it will still increase to a certain extent, which may have a certain impact on Sinorhizobium. We should verify the tolerance threshold of Sinorhizobium to toxins through more sophisticated models in the future, which will be our future work.

Conclusion

Our models serve our program and are closely associated with the wet lab group. We modeled the formation of nodules and compared them to grown nodules, finding some similarities. Besides, we confirmed the broad spectrum of our project through molecular docking modeling. Additionally, with DHA yield prediction modeling, we made predictions for the outcome of our project, and the yield was reasonable and credible. Finally, the final suicide pathway ensures the safety of our project. We have corresponding modeling for each circuit to guide and improve the entire project.

The detailed codes of our models are all stored in GitLab.

References

Poole, P. Symbiosis for rhizobia is not an easy ride. (2024) Nat Microbiol 9, 314–315.
Wolfram, S. Statistical Mechanics of Cellular Automata. (1983) Reviews of Modern Physics, 55(3), 601-644.
Thompson, R. L., & Goel, N. S. (1988). Movable Finite Automata (MFA) models for biological systems. I: Bacteriophage assembly and operation. Journal of theoretical biology, 131(3), 351–385.
Jones, S., & Thornton, J. M. (1996). Principles of protein-protein and protein-DNA recognition. Proceedings of the National Academy of Sciences, 93(1), 13-20.
Yan, Y., Zhang, D., Zhou, P., & Li, H. (2017). HDOCK: a web server for protein-protein and protein-DNA/RNA docking based on a hybrid strategy. Nucleic Acids Research, 45(W1), W365-W373.
Yan Y, Wen Z, Wang X, Huang S-Y. (2017) Addressing recent docking challenges: A hybrid strategy to integrate template-based and free protein-protein docking. Proteins;85:497-512.
Huang S-Y, Zou X. (2014) A knowledge-based scoring function for protein-RNA interactions derived from a statistical mechanics-based iterative method. Nucleic Acids Res.;42:e55.
Huang S-Y, Zou X. (2008) An iterative knowledge-based scoring function for protein-protein recognition. Proteins;72:557-579.
Remmert M, Biegert A, Hauser A, Soing J. (2011) HHblits: lightning-fast iterative protein sequence searching by HMM-HMM alignment. Nat Methods.;9:173-175.
W. R. Pearson and D. J. Lipman. Improved tools for biological sequence comparison. (1988) Proc Natl Acad Sci USA;85:2444-2448.
Sievers F, Wilm A, Dineen D, Gibson TJ, Karplus K, Li W, Lopez R, McWilliam H, Remmert M, Soding J, Thompson JD, Higgins DG. (2011) Fast, scalable generation of high-quality protein multiple sequence alignments using Clustal Omega. Mol Syst Biol.;7:539.
Larkin MA, Blackshields G, Brown NP, Chenna R, McGettigan PA, McWilliam H, Valentin F, Wallace IM, Wilm A, Lopez R, Thompson JD, Gibson TJ, Higgins DG. (2007) Clustal W and Clustal X version 2.0. Bioinformatics.;23:2947-8.
Marti-Renom MA, Stuart A, Fiser A, Sanchez R, Melo F, Sali A. (2000) Comparative protein structure modeling of genes and genomes. Annu Rev Biophys Biomol Struct;29:291-325.
Berman HM, Westbrook J, Feng Z, Gilliland G, Bhat TN, Weissig H, Shindyalov IN, Bourne PE. (2000) The Protein Data Bank. Nucleic Acids Res;28:235-242
Resendis-Antonio, O., Reed, J. L., Encarnación, S., Collado-Vides, J., & Palsson, B. Ø. (2007). Metabolic reconstruction and modeling of nitrogen fixation in Rhizobium etli. PLoS computational biology, 3(10), 1887–1895.
Monk, J. M., Koza, A., Campodonico, M. A., Machado, D., Seoane, J. M., Palsson, B. O., Herrgård, M. J., & Feist, A. M. (2016). Multi-omics Quantification of Species Variation of Escherichia coli Links Molecular Features with Strain Phenotypes. Cell systems, 3(3), 238–251.e12.
Syska, C., Kiers, A., Rancurel, C., Bailly-Bechet, M., Lipuma, J., Alloing, G., Garcia, I., & Dupont, L. (2024). VapC10 toxin of the legume symbiont Sinorhizobium meliloti targets tRNASer and controls intracellular lifestyle. The ISME journal, 18(1), wrae015.
Lipuma, J., Cinege, G., Bodogai, M., Oláh, B., Kiers, A., Endre, G., Dupont, L., & Dusha, I. (2014). A vapBC-type toxin-antitoxin module of Sinorhizobium meliloti influences symbiotic efficiency and nodule senescence of Medicago sativa. Environmental microbiology, 16(12), 3714–3729.
Bugara, B., Konieczny, P., Wolnicka-Glubisz, A., Eckhart, L., Fischer, H., Skalniak, L., Borowczyk-Michalowska, J., Drukala, J., & Jura, J. (2017). MCPIP1 contributes to the inflammatory response of UVB-treated keratinocytes. Journal of dermatological science, 87(1), 10–18.
Part:BBa K3512042 - parts.igem.org. [2024-09-25]. https://parts.igem.org/Part:BBa_K3512042.

Parameters	Description	Value	Unit
\(Ac_{in}\)	Input Concentration of Acetyl-CoA	\(1\)	\(\text{mol/L}\)
\(Ac_{out}\)	Output Concentration of Acetyl-CoA	\(1\)	\(\text{mol/L}\)
\(diff\)	Diffusion Rate	\(0.1\)	\(\text{mol/(L·h)}\)
\(v_{PHBmax}\)	Maximum PHB Production Rate	\(1\)	\(\text{mol/(L·h)}\)
\(K_{mPHB}\)	PHB Michaelis Constant	\(0.5\)	\(\text{mol/L}\)
\(v_{EPAmax}\)	Maximum EPA Production Rate	\(1\)	\(\text{mol/(L·h)}\)
\(K_{mEPA}\)	EPA Michaelis Constant	\(0.5\)	\(\text{mol/L}\)
\(c_{Ac0}\)	Initial Acetyl-CoA Concentration	\(0.1\)	\(\text{mol/L}\)
\(c_{PHB0}\)	Initial PHB Concentration	\(0.1\)	\(\text{mol/L}\)
\(c_{EPA0}\)	Initial EPA Concentration	\(0.1\)	\(\text{mol/L}\)
\(DHA_{c0}\)	Initial DHA Concentration	\(0.0\)	\(\text{mol/L}\)
\(DHA_{e0}\)	Initial DHA Secretion Concentration	\(0.0\)	\(\text{mol/L}\)
\(k\)	DHA Synthesis Rate Constant	\(0.5\)	\(\text{mol/(L·h)}\)
\(K_{DHA}\)	DHA Synthesis Rate Threshold	\(5.0\)	\(\text{mol/L}\)
\(S_{max}\)	Maximum DHA Secretion Rate	\(0.5\)	\(\text{mol/(L·h)}\)
\(K_{S}\)	DHA Secretion Michaelis Constant	\(2.0\)	\(\text{mol/L}\)

Parameters	Description	Value	Unit
\(v_{M_C}\)	\(M_C\) transcription rate	\(0.1333\)	\(s^{-1}\)
\(v_{M_B}\)	\(M_B\) transcription rate	\(0.1333\)	\(s^{-1}\)
\(v_{sg}\)	\(sg\) transcription rate	\(0.1333\)	\(s^{-1}\)
\(v_{M_K}\)	\(M_K\) transcription rate	\(0.1333\)	\(s^{-1}\)
\(\beta_{C}\)	\(M_C\) rate of translation	\(0.033\)	\(s^{-1}\)
\(\beta_{B}\)	\(M_B\) rate of translation	\(0.139\)	\(s^{-1}\)
\(\beta_{K}\)	\(M_K\) rate of translation	\(0.139\)	\(s^{-1}\)
\(\delta_{M_C}\)	\(M_C\) rate of decay	\(0.00203\)	\(s^{-1}\)
\(\delta_{M_B}\)	\(M_B\) rate of decay	\(0.00203\)	\(s^{-1}\)
\(\delta_{M_K}\)	\(M_K\) rate of decay	\(0.00203\)	\(s^{-1}\)
\(\delta_{sg}\)	\(sg\) rate of decay	\(0.00203\)	\(s^{-1}\)
\(d_C\)	\(C\) rate of decay	\(2.8881 \times 10^{-4}\)	\(s^{-1}\)
\(d_B\)	\(B\) rate of decay	\(0.00115524\)	\(s^{-1}\)
\(d_K\)	\(K\) rate of decay	\(0.00115524\)	\(s^{-1}\)
\(F_{KC}\)	Cas12k rate of inhibition	\(0.2\)	-
\(a_{CB}\)	CB rate of complex formation	\(2 \times 10^{-3}\)	\(M^{-1} s^{-1}\)
\(\theta_{sgC}\)	\(sg\) binding C-gene rate	\(7.14972 \times 10^{-3}\)	\(M^{-1} s^{-1}\)
\(K_D\)	\(glnK-N\) response index \(nifH-O_2\) response index	\(800 \times 10^{-6}\)	\(s^{-1}\)