Summary
Modeling is a fundamental part of the BioPank project since evaluating the paratransgenesis approach in the field or in the laboratory would require controlling and measuring different conditions. As detailed in this section, We built a kinetic model with parameter optimization, where an Ordinary DIfferential Equations (ODEs) framework coupled with sensitivity analysis evaluates the performance of our chassis under variable conditions in the midgut of the vector and of the engineered genetic circuit to eliminate leishmania. Also, we built the Carpincho Toehold Generator, specialized in the design of Toeholds for miRNAs detection and an essential tool to design our detection module throughout the project.
Kinetic Model
Motivation
The simulation exploits a preliminary screenshot of the circuit's behaviour to guide future experiments in the lab. In general, the persistence of the chassis and the efficiency of eliminating the parasite within the midgut are essential factors in paratransgenesis. For the success of our approach, we needed to understand how different conditions affect these events. Therefore, this kinetic model aimed to visualize leishmania death and chassis growth behavior over time. To optimize the system, we identified the parameters that control AMP production, chassis growth and parasite death in order to maximize i) the time needed to produce sufficient AMP to eliminate the parasite and ii) the chassis persistence within the vector.
Model Structure
This mechanistic model contains Detection and Elimination modules influenced by global transcriptional regulation based on chassis growth in the vector's midgut. The mass action kinetics describe the kinetic model applied to chassis growth, Toehold transcription, AMP translation, parasite elimination and degradation of molecules, resulting in seven entities (substrate [S], chassis [Bs], Toehold mRNA [mRNAtoehold], activated Toehold mRNA [mRNAon], proAMP [AMPoff] and AMP [AMPon] concentrations). The equations, involving 18 parameters (described in Table 1), describe the 15 reactions of the entire process. In the following representation, we visualize the relationships between parameters and variables
Figure 1: Representation of reactions and biological mechanisms in our system.
Based on this representation, we divided the system into four parts:
- Growth Model, which describes Bacillus subtilis [Bs] concentration as a function of substrate availability (Monod kinetic equation) and concentration decay as a function of re-sporulation probability, decay rate due to varied conditions; and cytotoxic feedback of antimicrobial peptides [AMP]:
$$\frac{d[Bs]}{dt} = V_{Bs,growth}-V_{Bs,death}$$
- Substrate Flux, which indicates the substrate input (vector digestion) as a function of cosine wave, chassis consumption and substrate removal [S] over time in the midgut of the vector:
$$\frac{d[S]}{dt} = V_{Bs,intake}-V_{S,rem}$$
- Antimicrobial Peptide Expression Model (color), which describes the constitutive transcription of Toehold mRNA [mRNAtoehold] as a function of chassis growth rate; the Toehold binding to miRNA trigger [mRNAon], considering the thermodynamic values obtained in the Carpincho tool (see below), enabling the proAMP translation [AMPoff]; and AMP acitvation [AMPon], through its diffusion to the extracellular medium and enzymatic digestion by trypsin at a constant rate. Also, we performed a rough approximation by multiplying the single cell AMP production per chassis concentration (cell/L) to measure the leishmanicidal activity. The AMPs and RNAs undergo first order degradation rates:
$$\frac{d[mRNA_{toehold}]}{dt} = V_{tx,yoehold}-(V_{mRNA,on}+V_{mRNA,deg})$$
$$\frac{d[mRNA_{on}]}{dt} = V_{mRNA,on}-V_{mRNAon,deg}$$
$$\frac{d[AMP_{off}]}{dt} = V_{tl,AMP}-(V_{act,AMP}+V_{AMP,deg})$$
$$\frac{d[AMP_{on}]}{dt} = V_{act,AMP}-V_{AMP,deg}$$
- Leishmania Death Model, which describes the influence of the concentration of activated antimicrobial peptides [AMPon] on leishmania mortality as a function of Hill's equation:
$$\frac{d[Leish]}{dt} = V_{Leish,growth}-V_{Leish,kill}$$
Major assumptions, before starting the simulation and analysis, we highlight the major assumptions considered in this model:
Substrate input has oscillatory behavior with the same frequency and amplitude;
The presence of the miRNA trigger is constant;
We only evaluate the chassis entry moment into the vector;
For more assumptions, check below our complete model description.
Complete Code
The Simulations and sensitivity analysis were performed using Colab Notebook for Python. The simulation graphs were plotted using the "plotly" package and the sensitive analysis graphs were plotted in RStudio. Codes and Images can be downloaded here.
Simulation Evaluation
To assess the feasibility of our project, we solved the single system of ordinary differential equations using the "odeint" function from the SciPy package for Python. Each reaction has been built upon parameter and initial condition regarding appropriate kinetic rate laws at steady state (Table 1). The simulation presented the results of the dynamics of our chassis concentrations, the substrate flux, the expression circuit, AMPs concentrations and parasite death.
The prediction is consistent with the dynamics of our system (Fig. 2,3 and 4) and the results were positive for eliminating leishmania. In Figure 2, as expected using the Monod equation, which considers the concentration of a nutrient-limited (in this case, glucose), the chassis colony had a growth directly related to the substrate input rate. On the other hand, it is observed that the transcription and translation had a disturbance of close to 200 minutes. This oscillation happens because we describe transcription as a function of global regulation, where growth rate is directly related to the constitutive promoter activity and dilution rates [9,15]. The drop of substrate availability decreases the growth rate and concentration of viable chassis, which induces a negative perturbation in the expression levels in Figures 3 and 4.
Figure 2. Dynamic of the concentration of our chassis (A) and substrate (B).
Figure 3. Dynamic of single-cell mRNA concentrations within 500 minutes (A) and 30 minutes of simulation (B) between activated (green line) and deactivated (blue line) Toehold mRNAs. Dashed lines: Time to Half maximal Toehold activation from transcripts.
Figure 4. A) Dynamic of colony proAMP production (purple line) and activation(blue line). B) Dynamic of Leishmania viable cells inside the midgut. Dashed lines: Time to reach half maximal inhibitory.
Even though this model integrated a strategy with the thermodynamic values obtained in the toehold design, its model neglects the uptake of miRNAs. In the model, the RNA uptake occurs at constant rates, which may have influenced the rapid activation shown in Figure 3, approximately 4 minutes for the activated toehold to reach half of its maximum concentration. Overall, these model simulations demonstrate the feasibility of our project to eliminate the parasite within the vector. The time needed to reach half-maximal inhibitory behavior was around 115 minutes (Fig. 4B) and keeping these positive results, we decided to go deeper and analyze the influence of each parameter on our system.
Analysis and Parameter Estimation
We could not measure all parameters experimentally, especially due to difficulties in lab access this year and the number of parameters required to evaluate a paratransgenesis system. As shown in Table 1, we used parameter values found in literature, which may not match our scenario since biological systems and experiments vary considerably between studies [17]. Furthermore, biological systems are subject to stochasticity and noise, which can make a parameter estimation difficult for nonlinear models [18]. We applied a Global Sensitivity Analysis (GSA) that seeks to determine the relative effect of the main reactions/parameters on the system-wide output dynamics through simultaneous perturbations on all inputs during several samplings.
Global Sensitivity Analysis
The lack of computational capacity and time available motivated our team to perform this analysis only in the last time-step of the simulation (last step = 500min). We set up the GSA with 520,000 parameter sets (from 10,000 samples for each variable and 25 parameters, sets = N(2D+2)), randomly sampled within 5% of their reference values (Table 1). As a result, we obtained the values of the first-, second and total-order Sobol Indices to analyze the significance of each parameter for AMP production, chassis growth and parasite death.
In Figure 5, the First order sensitivity indices has six values greater than 0.2 for each component analyzed. This observation demonstrates that some parameters alone have the potential to affect the system significantly. Here, we consider values close to or above 0.4, such as n (AMPs Hill Coefficient) and v (B. subtilis maximum growth rate).
On another hand, second order sensitivity indices show interactions between the parameters with a low level of effect on the model output, as shown in Figure 6. The interaction between v and lamb (single-cell half-life for B. subtilis persistence decay rate, Table 1) will be negligible due to the low index value.
Figure 6. Second Order sensitivity indices for the last simulation step (500min).
In Figure 7, Total-order sensitivity indices contain the same aspect as the first-order indices. However, the lamb index probably increases due to its interaction analyzed in Figure 6.
Total-Order Sobol Indices
Figure 7. Total Order sensitivity indices for the last simulation step (500min).
Overall, the model and analysis achieve important directions despite the system not showing robustness in specific parameters. First, we note that the amount of AMP is significantly affected by gene expression circuit parameters that describe the degradation rates of the molecule (i.e., DEGrna), transcription, and enzymatic trypsin digestion (i.e., dig_rat). It is noteworthy that these reactions have been simplified and most of the variables related to their mechanism have been neglected (see Complete Description).
Second, the concentrations of activated AMP, chassis and parasite are directly influenced through kinetic characteristics related to chassis growth and AMP cooperativity. These values are consistent with their interaction throughout the system, as the cooperativity of AMP (n) can affect the growth of our chassis (remember that AMP can be toxic for both parasite and chassis) and, consequently, the production of AMP. Also, the maximum growth (v) rate can interfere with the transcription rate and directly affect the concentration of Leishmania (remember that we describe our transcription rate with global regulation).
To assess their impact on the time response to reach 50% inhibition of Leishmania viable cells, we simulated a wide range of AMP cooperativity, n (1 to 4), and a wide range of maximum growth rates for the chassis, v (0 to 0.064 min-1). In figure 8, we see that the v parameter strongly affects the response time for half elimination. Furthermore, it is important to consider that we only assume the analysis starting from the chassis input in the vector, that is, considering the proAMPs system, we could have previously high AMP levels.
Figure 8. Half-Elimination Response time as function of two significantly parameters: n (AMP's Hill coefficient) and v (maximum growth rate of Bacillus subtilis)
Finally, despite showing a non-robust behavior in the variation of these parameters, efforts to optimize the system and the quick response can be focused on the engineering of AMPs and on the growth capacity of our chassis. We can evaluate different AMPs with different cooperativity mechanisms and also optimize substrate uptake by the chassis. Furthermore, we note that the chassis decay rate within the vector (lamb) must be carefully evaluated in future experiments (see Protection Module). Its interaction with chassis growth properties significantly affects the output system dynamics.
Conclusions
Our model and analysis show promising system efficiency results and helped obtain screenshots to guide future efforts. The robustness of the system must be evaluated to tune the most significant parameters. The way we built the model helped in these future steps, as we focused on measurable parameters for wet lab analysis conditions that we have in our laboratory. Hill's coefficient (n) and maximum growth rate (v) appeared as the most relevant parameters to be tuned, leading to experiments trying to assess Hill coefficient for our AMPs candidates during this season (see Experiments). Finally, more robust analyses can increase system understanding. For instance, some negative values appear in the GSA and we evaluate the output of a time-step. The Sobol method considers unusual negative values, and the GSA must be evaluated along all time-steps due to the simulation non-linearity. More extensive sampled analyses can solve these factors.
Carpincho Toehold Designer
The design of toehold switches requires an in silico knowledge of the thermodynamic characteristics that can influence their application. Our team built the Carpincho Toehold Designer, a slow but efficient way to design our Toeholds for detecting miR-283 (from L. longipalpis) obtained from midgut by Bacillus subtilis. The Carpincho Tool strictly follows the Type II model described by Green et al. (2020) and bases its operation on NUPACK Design and Analysis algorithms [20].
Carpincho is a Toeholds generator that requires low computational knowledge and requires few inputs to generate the libraries. We calculate the Minimum Free Energy (MFE) variation between the binding and unbound state through the energy difference between the structures formed by Trigger and Toehold binding with the MFE of conformation of Trigger and Toehold alone. Then, we can select the best candidates according to the relationship:
$$\Delta G_{DeltaMFE}<\Delta G_{Switch}+\Delta G_{mRNA Trigger}$$
In addition, other factors can be considered, such as ensemble_defects and MFE RBS_Linker, parameters that can optimize Toehold activation and decrease the leaky expression. The code and the generated results (conformation figures and thermodynamic data) can be consulted in our repository. Despite having some limitations, Carpincho can be optimized to expand its design potential to other miRNAs in different organisms.
References
[1]. Ugalde-Salas, P. et al. Insights from Microbial Transition State Theory on Monod's Affinity Constant. Nature Scientific Reports, 2020.
[2]. Heerman, M. et al. Bacterial Infection and Immune Responses in Lutzomyia longipalpis Sand Fly Larvae Midgut. PLOS Neglected Tropical Diseases, 2015.
[3]. Gauvry, E. et al. Differentiation of Vegetative Cells into Spores: a Kinetic Model Applied to Bacillus subtilis. Applied and Environmental Microbiology, 2019
[4]. Wang, S. et al. Driving mosquito refractoriness to Plasmodium falciparum with engineered symbiotic bacteria. Science, 2017.
[5]. BioNumbers ID 115203
[6]. Wang, X. Modeling of the Bacillus subtilis Bacterial Biofilm Growing on an Agar Substrate. Hindawi: Computational and Mathematical Methods in Medicine, 2015.
[7]. Overkamp, W. et al. Physiological and cell morphology adaptation of Bacillus subtilis at near-zero specific growth rates: a transcriptome analysis. Environmental Microbiology, 2014.
[8]. Gauvry, E. et al. Differentiation of Vegetative Cells into Spores: a Kinetic Model Applied to Bacillus subtilis. Applied and Environmental Microbiology, 2019
[9]. Gerosa, L. 2013. Dissecting specific and global transcriptional regulation of bacterial gene expression. Molecular Systems Biology, 2013.
[10]. Brand, G. D. et al. The Skin Secretion of the Amphibian Phyllomedusa nordestina: A source of Antimicrobial and Antiprotozoal Peptides. Molecules, 2013.
[11]. Brand, G. D. et al. Novel dermaseptins from Phyllomedusa hypochondialis (Amphibia). Biochemical and Biophysical Research Communications, 2006.
[12]. Pérez-Cordero, J.J. et al. (2011) Leishmanicidal activity of synthetic antimicrobial peptides in an infection model with human dendritic cells. Peptides, 32(4), p. 683-690
[13]. Serafim, T. D. et al. Sequential blood meals promote Leishmania replication and reverse metacyclogenesis augmenting vector infectivity. Nature Microbiology, 2018
[14]. Radeck, J. et al. The Bacillus Biobrick Box: generation and evaluation of essential genetic building blocks for standardized work with Bacillus subtilis. Journal of Biological Engineering. 2013.
[15]. Hintsche, M.;Klumpp,S. (2013) Dilution and theorical description of growth-rate dependent gene expression. Journal of Biological Engineering, 7(22).
[16]. Senoussi, A. Quantitative characterization of translational riboregulator using an in vitro transcription-translation system. ACS Synthetic Biology, 2018.
[17]. Kent, E. et al.(2013) What Can We Learn from Global Sensitivity Analysis of Biochemical Systems?. PLoS ONE 8(11):e79244.
[18]. Quang, M. N. et al. (2019) Global Sensitivity Analysis of Metabolic Models for Phosphorus Accumulating Organisms in Enhanced Biological Phosphorus Removal. Frontiers in Bioengineering and Biotechnology. 237(7).
[19]. Green, A.A. et al. (2020). Complex cellular logic computation using ribocomputing devices. Nature, 548, p. 117-121.
[20]. Green, A.A. et al. (2020). Complex cellular logic computation using ribocomputing devices. Nature, 548, p. 117-121.