Why Agent-Based Models?
Predicting how genetically engineered P. fluorescens behaves in soil is complicated due to intricate microbial interactions and fluctuating environments. This is not an easy problem to solve for deterministic models, as the actions of a single bacteria/colony might lead to unknown collective behaviors. With several simple rules for bacteria to follow, agent-based models can explore complex biological systems on multiple levels [1], including entity behaviors, microbial interactions and environmental effects (Figure 1).
Figure 1. General scheme of agent-based modelling. Inputs from internal or external factors and stochastic choices made by agents eventually form the outputs of the model.
Aim of the Model
To predict more precisely the behaviors of our genetically engineered P. fluorescens in the soil, we have built two sub-models: a Microbial Community model to simulate how our bacteria move and colonise on the root; and a Cell Level model to explore how the protein secretion system would work in cells with genetic circuits. Addtionally, we intertwined two sub-models to explore cell-level effects on antiflorigen protein secretion across populations, in pursuit of the ideal product application. Eventually, we are trying to give suggestions for both experimental design and product development. Here, we are using NetLogo [2] as an agent-based modeling platform to build our model and conduct simulations; this also allows us to even track the spatial movements of the entities during the simulations.
What We Learned?
From our agent-based models, we got new insights into the behaviours of our bacteria in the soil. This helps us in both experiments and product development. Moreover, we also integrated the results of our agent-based models to the Applied Design, giving suggestions for the product design and precisely application methods. Furthermore, through the analysis of large amounts of data from the agent-based models, we successfully captured a series of real-life biological patterns. This proves that agent-based modelling frameworks can be potentially powerful tools for the complex biological systems.
Microbial Community Model
In the Microbial Community model, we explore the bacteria behaviours at population level in the soil. The agents in this sub-model are natural P. fluorescens and our modified P. fluorescens. Those agents will follow several simple rules (see below) to act in a 20 x 20 cm closed field in the soil. The output of this sub-model is colonisation efficiency of our modified P. fluorescens on roots (i.e., the percentage of population of our bacteria on the root).
You can also enjoy the non-embedded Microbial Community model here in a separate page.
For smaller screens this content is not embedded. Instead go and enjoy the microbial community model here.Model Validation
Figure 2. Spatial distribution of colonisation of modified P. fluorescens on roots. (a) Colonisation distribution of our bacteria on roots with smaller root size. (b) Colonisation distribution of our bacteria on roots with bigger root size. (c) Part of experimental results from previous research on root colonisation of P. fluorescens [11]. (d) Our experimental measurement of spatial distribution of colonisation of modified P. fluorescens on the root of Arabidopsis thaliana (Discover more in our WetLab page).
We modelled the spatial distribution of the modified P. fluorescens on the root by adding the population of our bacteria for each root region from 30 simulations; the darker colours are where more colonies are likely to be found (Figure 2a and 2b). Given the comparable root lengths between wheat and cherry trees, we utilized wheat data for a comprehensive comparison (Figure 2c). We can observe that our bacteria will colonise mostly in the first 6 cm underground on root, and the colonisations are showing a decreasing gradient from top layer to deeper layer of the roots in the soil (Figure 2a and 2b), as shown in the wheat root (Figure 2c). Moreover, we experimentally measured the colonisation patterns of our bacteria on the root of Arabidopsis thaliana. Our model's findings corresponded with the observed data from the wet lab (Figure 2c). This means the model can capture the essential aspects of the real-world system.
Colonisation Dynamics
Figure 3. Example plots of colonisation dynamics of P. fluorescens on the root. 4 graphs show 4 different parameter conditions to the system. (a) Water-content-of-soil: 14%; root-cc: 200000; competition-death-rate: 0.014 /h; max-cell-growth-rate: 0.4 /h. (b) Water-content-of-soil: 30%; root-cc: 200000; competition-death-rate: 0.014 /h; max-cell-growth-rate: 0.5 /h. (c) Water-content-of-soil: 22%; root-cc: 400000; competition-death-rate: 0.25 /h; max-cell-growth-rate: 0.5 /h. (d) Water-content-of-soil: 22%; root-cc: 500000; competition-death-rate: 0.014 /h; max-cell-growth-rate: 0.5 /h. The blue line shows the population dynamics of the natural P. fluorescens with error bars. The orange line shows the population dynamics of the modified P. fluorescens with error bars. The red line is the standard line of 50% of colonisation efficiency for comparing the dynamic patterns.
We conducted a parameter sweep, exploring combinations of parameter values through an extensive series of simulations. Here we show four example plots from different parameter conditions with 30 repetitions for each condition, which leads to different colonisation dynamic patterns (Figure 3). Within the simulation time, the population of natural P. fluorescens and modified P. fluorescens will reach steady states, or constant levels, after around four to nine days. This means if we want to apply our product to the orchard, it would be better to apply it at least ten days earlier to make sure our product works for the plants (see Model Suggestions below).
Predictability of the Microbial Community Model
Figure 4. Random Forest regression for Microbial Community model. The blue dots are the output values (i.e., colonisation efficiency) from both trained data and predicted data. The red line is the regression line of the data points.
Subsequently, we used random forest regression to test the predictability of the Microbial Community model. Despite being a stochastic model, the mean squared error is low and the R-squared score is over 0.99 (Figure 4), which means the behavior of our bacteria would be relatively controlable at the population level in the soil.
Cell Level Model
Figure 5. General scheme of the genetic circuits.
The biosensor system comprises two genetic circuits: one activated by a specific root exudate – salicylic acid (click here to find out the details of screening of root exudates in weblab page) and the other by N-Acyl-l-homoserine lactone (AHL) (Figure 5). AHL is continually produced but functional to the cells only above a certain concentration threshold. Consequently, the target protein is synthesized only in the presence of salicylic acid and a sufficiently bacterial population.
In the Cell Level model, we simplify several genetic circuits into only one, and the production of trigger RNAs and AHL is expressed in simpler ways. In this sub-model, we assume that AHL concentration is always high enough to activate the protein secretion. The agents of this sub-model are the molecules in the cell (e.g., RNA polymerase, proteins, etc.). The setting of the strength of promoter, ribosome binding site, and protein degradation rate are mostly set to default following an example model in NetLogo [12]. This sub-model will follow several rules (see below) to perform transcription, translation, etc., exploring the protein secretion of cells with artificial genetic circuits.
You can also enjoy the non-embedded Cell Level model here in a separate page.
For smaller screens this content is not embedded. Instead go and enjoy the cell level model here.Model Validation
Figure 6. Comparison of Cell Level model’s results and lab data. “R + AHL” represents the condition that both root exudate salicylic acid and AHL occur in the system. “Leaky Expression” represents the condition that the target gene is not activated in the system. The blue and red lines are experimental measured data. The orange and black lines are model simulations with error bars.
We compared our model with our own generated lab data; the model’s results were collected after 30 times repetitions of simulation (orange and black regions) and shows that the results from the Cell Level model fits the data generated from lab (blue and red lines) (Figure 6). This means the Cell Level model can explain the protein secretion system in the cells.
Predictability of the Cell Level Model
Figure 7. Random Forest regression for Cell Level model. The blue dots are the output values from both trained data and predicted data. The red line is the regression line of the data points.
Unlike the Microbial Communities model, the Cell Level model shows more variability due to random movements of molecules in the cells, making the protein secretion more unpredictable based on the parameters. This is reflected by the larger mean squared error and lower R-squared value (Figure 7). We hypothesise that this is because the dynamics of protein production could be different even though the basic trend of protein secretion is like the real data.
Model Interaction
We used a built-in tool LevelSpace [13] in NetLogo to conduct model interactions. The output of the interactive model is the total protein secretion level in the system.
Rule and kinetic model for the interactive model:
Rule:
The cell level model will be activated only when the AHL concentration for colonies are above the activation threshold [14].
Kinetic model:
AHL production model follows the equation from previous research [15]:
$katex$\dfrac{dA}{\phantom{.}dt}= - \gamma A + \dfrac{\alpha + \beta A^{n}}{A_{threshold}+A^{n}}\cdot{} N(t)$katex$
where α is constitutive production rate of AHL. β is induced production rate of AHL. n is the degree of polymerization. Athreshold represents the threshold of AHL between low and increased activity. γ is the degradation rate of AHL. N(t) is the population of the bacteria, changing with time.
Watch the video to see how the interactive model works:
Video 1: An example simulation of the interactive model. The built-in tool - LevelSpace does not support to run in the browers. So, here we show a video of the model interaction.
If the video is not loading, watch it here.
Predictability of the Interactive Model
Figure 8. Random Forest regression of the interactive model. The blue dots are the output values from both trained data and predicted data. The red line is the regression line of the data points.
We integrated the protein secretion level of individual cells with the population of colonised genetically modified P. fluorescens, thereby encompassing both colonisation efficiency and protein secretion magnitude. Perturbations from the Cell Level model affects the predictability of the interactive model, which has R-squared score of 0.46 (Figure 8) compared to 0.99 (Figure 4) before. Most of the errors likely come from stochasticity of the sub-models.
Sensitivity Analysis
Sensitivity analysis is crucial as it reveals the impact of parameter variations on mathematical models. It identifies key parameters, guiding experimental focus and resource allocation. By pinpointing influential factors, it enhances model accuracy, aiding in decision-making and advancing scientific understanding.
SHapley Additive exPlanations (SHAP)
SHAP is a data analysis method that explains how each feature impacts the model's output. It compares the model’s output for all possible variable combinations both with and without one of the variables [16]. Here, we use SHAP values to conduct sensitivity analysis at both population and cell levels.
Figure 9. Sensitivity analysis for the agent-based model. (a) and (c) demonstrate the SHAP analysis of variables in population level and cell level, respectively. (b) and (d) demonstrate the feature importance of each variable by calculating the mean absolute SHAP values.
Root size (“root-cc”) is the most sensitive parameter in the interactive model. This is likely because bigger roots have larger carrying capacity for bacteria. The second most sensitive parameter is competition death rate, which significantly influences total protein secretion. Therefore, enhancing bacterial competition ability isn't necessary, allowing us to control costs effectively. We also discovered that properly increasing maximum cell growth rate and water content of soil would also increase the total protein secretion level (Figure 9a and 9b). Those findings provide deeper insights into product development and application strategies as will be discussed in model suggestions.
The protein degradation rate is the most sensitive parameter in the cell level model, which means that the stability of our target protein is vital for the protein secretion (Figure 9c and 9d). Our results suggest that it is advantageous to perform laboratory evolution for the enzyme of our bacteria to increase the stability of the desired protein in the future.
Model Suggestions
Using a reported standard curve of relationship between sGFP fluorescence and sGFP protein concentration [17], we replot a calibration curve to estimate the range of total protein concentration of the final system (Figure 10).
Figure 10. sGFP standard curve. The blue dots are the values of sGFP fluorescence related to the sGFP concentrations [17]. The orange line is linear regression curve of these data points.
Based on the data generated from different combinations of parameter values (Figure 11), the result shows that the range of total protein concentration after 2 weeks’ simulation can be at least 33.3 mg/ml, and up to 6510 mg/ml in optimized condition.
Figure 11. Visualization of all data points of the interactive model’s output. The green dot is the model output with the highest protein secretion level. The red dot is the model output with the lowest protein secretion level. The orange and purple dots are the model outputs with fixed condition, except 14% water content of the soil and 30% water content of the soil respectively. The blue dots are other data points of the model outputs.
The introduction of florigen in Arabidopsis necessitates concentrations ranging from tens to hundreds of micromoles [18]. In contrast, the maintenance of flowering balance typically requires relatively small amounts of antiflorigen factors [19]. So, we provide far more than enough antiflorigen protein to restore the flowering of plants. Therefore, we can have a lower proportion of our bacteria in the product to control the cost.
The outcome demonstrates our product's effectiveness in controlling plant flowering. If necessary, enhancing total protein secretion can be achieved by improving bacterial growth or protein stability. Additionally, maintaining well-hydrated soil is advised for better product performance in agricultural applications. Although our model says we need to know temperature four days in advance, given we are overproducing antiflorigen protein, we believe that earlier application would also be beneficial.
