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
- imagem -
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 Model Description
Our complete ODEs system's description and assumptions can be checked in the next file. Also, 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 can be downloaded [here].
Model parameters were collected from publications, mainly because it would require several experiments to quantify the values in detail. Still, some experimental results were used to determine essential parameters such as inhibitory concentrations of AMPs against Leishmania infantum, check our Results. The table footnotes describe some estimates or rescaling calculations.
| Name | Value | Unit | Source | Description | Gsa Notation |
|---|---|---|---|---|---|
| Growth Model | |||||
| $$\mu$$ | Substrate-dependent | h^-1 | [1] | B. subtilis growth rate. | - |
| $$t_{\frac{1}{2}}$$ from b | 1.5(90) | h^-1 (min) | [2] | single-cell half-life for B. subtilis persistence decay rate | lam |
a) Mean values of B. subtilis persistence within 24h in L. longipalpis larvae gut;
b) 10^2 Cells/L = 10^5 CFU/mL from [4]
c) Bacteria unit per substrate concentration (rescaled)
d) This value remains for each substrate uptake over time (See Substrate Flux).
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.
asd²