Toxin-Antitoxin Model

Overview

In this model, we attempt to use experimental data to determine the relations between the toxin (Ccdb) and antitoxin (Ccda) in our root-localization kill switch system.

Our root-localization kill switch system, based on malic acid in the root exudate, expresses anti-toxin under the influence of malic acid (malate). Under ideal conditions, the antitoxin production rate should be able to outweigh the production of the toxin (given enough malic acid); however, if the bacteria travels too far away from the roots, the concentration of malic acid decreases, causing the antitoxin's production to decrease as well. The decrease in the antitoxin production would kill the bacteria, preventing it from spreading away from the banana plant.

Our purpose in modeling this system is to determine the survivability of the bacteria given the concentration of the malic acid and by extension, how our solution should be implemented.

The complete construction and design of the root-localization kill switch can be found here: Wetlab Page of Experiments. Essentially, the malate binds with a constitutively expressed activator (mleR) and increases the translation of CcdA. CcdB is also constitutively expressed, maintained at a low concentration that can be counteracted by CcdA, but only when malic acid increases CcdA expression.

CcdA Reactions

CcdA Reactions

To create the mathematical equations to model the system, we first organized it into three separate reactions, labeled as 1, 2, and 3.

\begin{equation} \begin{align*} \ce {[\text{MA}] +[\text{mleR}] <=>[k_{on_1}][k_{off_1}] [MM]} \end{align*} \end{equation}

\begin{equation} \ce {[\text{pmle}] + [\text{MM}] <=>[k_{on}_2][k_{off}_2]} [\text{pmle} \times \text{MM}] \end{equation}

\begin{CD} \tag{3} [\text{pmle} \times \text{MM}] @>k_\text{mccda}>> [\text{mRNA}] @>k_{ccda}>> [\text{CcdA}] \\ @. @Vd_\text{mrna}VV @Vd_\text{ccda}VV \\ @. \phi @. \phi \\ \end{CD}

Reaction 1 describes the relationship between malic acid (ligand) - [MA] and [mleR] (receptor), which combine to create a malic acid-mleR complex (activator), simply named [MM].

Reaction 2 describes the relationship between [MM] and [pmle] (mleR regulatory binding sites), creating [pmle x MM] (activated promoter).

Reaction 3 describes how [pmle x MM] determines the transcription of mRNA, which ultimately leads to the expression of CcdA.

From the reactions described above, we can create differential equations to model the system. From these equations, we can gain more insight into the dynamics of our system.

Reaction 1

Before we adapted reaction 1, we defined each of the terms in the reaction. $[\text{MA}]$ is the concentration of malic acid, $[\text{mleR}]$ is the concentration of mleR, $[\text{MM}]$ is the concentration of the malic acid-mleR complex, $\text k_{\text {on}_1}$is the association rate, and $\text k_{\text {off}_1}$ is the dissociation rate.

\begin{equation} \tag{1} \begin{align*} \ce {[\text{MA}] +[\text{mleR}] <=>[k_{on_1}][k_{off_1}] [MM]} \end{align*} \end{equation}

To adapt it into a mathematical equation (4), we utilized the Law of Mass Action, which describes the reaction rate as proportional to the concentration of each reactant:

\begin{equation} \tag{4} {d[\text{MM}] \over dt} = k_{\text{on}}[\text{MA}] - k_{\text{off}}[\text{MM}] \end{equation}

This gives us an ordinary differential equation relating to the rate at which $[\text{MM}]$ changes due to $[\text{MA}]$. Here, $\text k_{\text {on}_1}$ and $\text k_{\text {off}_1}$are constants that determine the “speed” of the equation. $\text k_{\text {on}_1}[\text{MA}]$ determines how much $[\text{MM}]$ is created at a given moment, while $\text k_{\text {off}_1}[\text{MM}]$ determines how $[\text{MM}]$ is reverted (dissociated) at a given moment.

Here, we assume that the $[\text{mleR}]$ is considered a constant and sufficient enough as to not affect translation in any meaningful way, so it is neglected in this equation.

Reaction 2

Like reaction 1, each term of reaction 2 is defined. $[\text{pmle}]$ is the concentration of pmle, $[\text{pmle} \times \text{MM}]$ is the concentration of the malic acid-mleR complex, $\text k_{\text {on}_2}$is the association rate of, and $\text k_{\text {off}_2}$ is the dissociation rate.

\begin{equation} \tag{2} \ce {[\text{pmle}] + [\text{MM}] <=>[k_{on}_2][k_{off}_2]} [\text{pmle} \times \text{MM}] \end{equation}

Using the Law of Mass Action, we derived reaction 2’s equation (5).

\begin{equation} \tag{5} {d[\text{pmle} \times {\text{MM}}] \over dt} = k_{\text{on}_2}[\text{pmle}][\text{MM}]- k_{\text{off}_2}[\text{pmle} \times \text{MM}] \end{equation}

This equation relates the rate of change of $[\text{pmle} \times \text{MM}]$ to both $[\text{pmle}]$ and $[\text{MM}]$.

Reaction 3

With the last reaction (3), we can split it into two smaller reactions (6 and 7) to simplify the creation of the mathematical model.

\begin{CD} \tag{3} [\text{pmle} \times \text{MM}] @>k_\text{mccda}>> [\text{mRNA}] @>k_{ccda}>> [\text{CcdA}] \\ @. @Vd_\text{mccda}VV @Vd_\text{ccda}VV \\ @. \phi @. \phi \\ \end{CD}

\begin{CD} \tag{6} [\text{pmle} \times \text{MM}] @>k_\text{mccda}>> [\text{mRNA}] @>d_{mccda}>> \phi \end{CD}

\begin{CD} \tag{7} [\text{mRNA}] @>k_\text{ccda}>> [\text{CcdA}] @>d_{ccda}>> \phi \end{CD}

For the transcription reaction (6), the terms are new defined as such: $[\text{mRNA}]$ is the concentration of CcdA’s mRNA, $\phi$ is degradation, $\text k_{\text {mccda}}$ is the transcription rate, and $\text d_{\text {mccda}}$ is the degradation rate of the mRNA.
Once again, using the Law of Mass Action, we can create a mathematical equation (8).

\begin{equation} \tag{8} {d[\text{mRNA}] \over dt}= k_{\text{mccda}}[\text{pmle}\times \text{MM}] - d_{\text{mccda}}[\text{mRNA}] \end{equation}

Here, the equation relates the rate of change of $[\text{mRNA}]$ to the $[\text{pmle} \times \text{MM}]$.
With these four equations as our building blocks, we created our final equation which will determine the concentration of CcdA given the concentration of malic acid.

QSSA Reduction

Due to the complexity of our model, we decided to simplify it by using Quasi Steady State Approximation (QSSA).The Quasi Steady State Approximation can be applied because the time scale of transcription (i.e. minutes) is far greater than that of translation (i.e. hours), so the rate-determining step of the reaction must be translation, the slower of the two. Thus, it can be simply assumed that malic acid rapidly binds with mleR and the resulting MM complex rapidly binds to the promoter and that the reaction reaches equilibrium in negligible time.

Reaction 1

We can apply QSSA assumptions to the equations we created above by substituting the rate of change to 0. Applying QSSA onto equation 4, we can solve for $[\text{MM}] $.

\begin{equation*} {d[\text{MM}] \over dt} = 0 \end{equation*}

\begin{equation*} 0 = k_{\text{on}_1} [\text{MA}]- k_{\text{off}_1}[\text{MM}] \end{equation*}

\begin{equation*} [\text{MM}] = {k_{\text{on}_1}[\text{MA}] \over k_{\text{off}_1}} \end{equation*}

Here, we simplified $\text k_{\text {on}_1}$and $\text k_{\text {off}_1}$by combining them into one variable, known as $\text K_{\text {d}_1}$, also known as the dissociation constant. In this case, $\text K_{\text {d}_1} = {k_{\text{off}_1} \over k_{on_1}}$.

\begin{equation} \tag{10} [\text{MM}] = {k_{\text{on}_1}[\text{MA}] \over k_{\text{off}_1}}= {[\text{MA}] \over K_{d1}} \end{equation}

Reaction 2

Moving on to the second equation (5), we manipulated the variables to solve for $[\text{pmle} \times \text{MM}]$.

\begin{equation} \tag{5} {d[\text{pmle} \times {\text{MM}}] \over dt} = k_{\text{on}_2}[\text{pmle}][\text{MM}]- k_{\text{off}_2}[\text{pmle} \times \text{MM}] \end{equation}

First, we wanted to remove $[\text{pmle}]$ from the equation to reduce the variables considered in the model. According to the equation shown below, we can do this since $[\text{pmle}]$ is the concentration of the free promoters.

\begin{equation*} \text{Free Promoter = Total Promoter - Bound Promoter} \end{equation*}

Utilizing this relationship, we can define Free Promoter as $[\text{pmle}]$, Total Promoter as Cn (the copy number of plasmids) and Bound Promoter as $[\text{pmle} \times \text{MM}]$.

\begin{equation} \tag{11} \text{[pmle] =} \:C_n - [\text{pmle} \times {\text{MM}}] \end{equation}

Substituting the equation (11) back into the ordinary differential equation (5), we get this:

\begin{equation*} {d[\text{pmle} \times {\text{MM}}] \over dt} = k_{\text{on}_2} (C_n-[\text{pmle} \times \text{MM}])[\text{MM}] - k_{\text{off}_2}[\text{pmle} \times \text{MM}] \end{equation*}

Doing another QSSA, we can assume the rate of change of the reaction reaches equilibrium, so ${d[\text{pmle} \times {\text{MM}}] \over dt} = 0$.

\begin{equation*} 0 = k_{\text{on}_2}(C_n-[\text{pmle} \times \text{MM}])[\text{MM}] - k_{\text{off}_2}[\text{pmle} \times \text{MM}] \end{equation*}

Now solving for $[\text{pmle} \times \text{MM}]$, we get equation 12:

\begin{equation} \tag{12} [\text{pmle}\times \text{MM}] = {[\text{MM}] \times C_n \over [\text{MM}] + K_{\text{d}_2}}, \text{where} \space K_{\text{d}_2} = {k_{\text{off}_2} \over k_{\text{on}_2}} \end{equation}

Reaction 3

Here with equation 8, we can simply do a QSSA to solve for $[\text{mRNA}]$ by assuming ${d[\text{mRNA}] \over dt} = 0$.

\begin{equation} \tag{8} {d[\text{mRNA}] \over dt}= k_{\text{mccda}}[\text{pmle}\times \text{MM}] - d_{\text{mccda}}[\text{mRNA}] \end{equation}

\begin{equation*} 0= k_{\text{mccda}}[\text{pmle}\times \text{MM}] - d_{\text{mccda}}[\text{mRNA}] \end{equation*}

Solving for $[\text{mRNA}]$, we get:

\begin{equation} \tag{13} [\text{mRNA}] = {k_{\text{mccda}}[\text{pmle} \times {\text{MM}}] \over d_{\text{mccda}}} \end{equation}

CcdA Final Equation

Lastly, with equation 9, we applied the rest of the equations that were simplified from their ODE forms.

\begin{equation*} \tag{9} {d[\text{CcdA}] \over dt}= k_{\text{ccda}}[\text{mRNA}]- d_{\text{ccda}}[\text{CcdA}] \end{equation*}

First, we substituted equation 10 into equation 12, which gives us this:

\begin{equation*} \tag{14} [\text{pmle}\times \text{MM}] = {[\text{MA}] \times C_n \over K_{\text{d}_1}\times K_{\text{d}_2}+[\text{MA}]} \end{equation*}

Then, we substituted equation 14 into equation 13, giving us this:

\begin{equation*} \tag{15} [\text{mRNA}] = {k_{\text{mccda}} \times C_n \over d_{\text{mccda}}} \times \Bigg({[\text{MA}] \over K_{\text{d}_1}\times K_{\text{d}_2}+[\text{MA}]}\Biggl) \end{equation*}

Now, substituting equation 15 into equation 9, we determined our final ODE that relates the concentration of malic acid to the concentration of CcdA.

\begin{equation*} {d[\text{CcdA}] \over dt} = {k_{\text{mccda}}\times {k_{\text{ccda}}} \times C_n \over d_{\text{mccda}}} \times \Biggl({[\text{MA}] \over K_{\text{d}_1}\times K_{\text{d}_2}+[\text{MA}]}\Biggl)- d_{\text{ccda}}[\text{CcdA}] \end{equation*}

\begin{equation} \tag{16} \end{equation}

Finally, we apply QSSA for the final time on equation 16 and consider the basal expression of CcdA ($\beta_0$). The basal expression is the base or leaky expression of a protein, essentially just the expression of CcdA without any malic acid in our case. We can solve for $[\text{CcdA}]$, giving us this:

\begin{equation*} \tag{17} [\text{CcdA}]= {k_{\text{mccda}}\times {k_{\text{ccda}}} \times C_n \over d_{\text{mccda}}\times d_{\text{ccda}}} \times \Biggl({[\text{MA}] \over K_{\text{d}_1}\times K_{\text{d}_2}+[\text{MA}]}\Biggl) +\beta_0 \end{equation*}

This is the equation on which our model is based. Noticeably, it is a Hill function with a Hill coefficient of 1, as malic acid and mleR bind in a 1:1 ratio, an independent binding.

CcdB Reactions

In addition to mathematically modeling CcdA, we also modeled the expression of CcdB. Since CcdB is produced constitutively, the reactions are much simpler than those of CcdA.

\begin{CD} [\text{pLac}] @>k_\text{mccdb}>> [\text{mRNA}] @>k_{ccdb}>> [\text{CcdB}] \\ @. @Vd_\text{mccdb}VV @Vd_\text{ccdb}VV \\ @. \phi @. \phi \\ \end{CD}

The reaction presented shows the gene expression of CcdB describing how [pLac] determines the transcription of mRNA, which ultimately leads to the expression of CcdB. Like before, using the Law of Mass Action, we can create two equations from this.

\begin{equation*} {d[\text{mRNA}] \over dt}= k_{\text{mccdb}}[\text{pLac}]- d_{\text{mccdb}}[\text{mRNA}] \end{equation*}

\begin{equation*} {d[\text{CcdB}] \over dt}= k_{\text{ccdb}}[\text{mRNA}]- d_{\text{ccdb}}[\text{CcdB}] \end{equation*}

We can then apply QSSA to reduce them from ODEs and substitute them with each other, arriving at this:

\begin{equation*} [\text{CcdB}]= {k_{\text{mccdb}}\times {k_{\text{ccdb}}} \times C_n \over d_{\text{mccdb}}\times d_{\text{ccdb}}} \end{equation*}

Assumptions

1. In our model, the concentration of mleR, the activator, was assumed to be fixed. The rate of mleR change is relatively negligible to total concentration. In our mathematical equation in equation 4, mleR concentration was assumed to be constant instead of a variable, only accounting for the change of malic acid concentration. In the lab, there are sufficient concentrations of mleR to binding with MA, making MA the rate-determining step. Therefore, the concentration of mleR is considered a constant and thus neglected in the mathematical equations. Small changes in the concentration of mleR will not meaningfully affect the rate of reaction.

2. We assumed that the malic acid is evenly distributed both within the cell and outside the cell. Since malic acid can diffuse through the cell wall and membrane, keeping the concentration balanced, we neglected a model describing the malic acid intake and usage and simply assumed malic acid concentration inside will equal the outside concentration. In addition, the malic acid is assumed to be released at a constant rate by the root.

3. We also assumed all proteins are distributed evenly in the cell so as to be unperturbed by binary fission.

4. Lastly, we assumed no outside factor impacting our bacteria, such as root turnover, weather, and other malic-acid consuming bacteria.

Finding Parameters

With both equations, we attempted to find each parameter of the model. We had very little success in finding the parameters in literature. As a result, we used parameters from previous iGEM projects for an estimate of most of the parameters. For the other ones, we estimated using literature estimates and protein modeling.

For estimated parameters, we utilized information from literature. For $\beta_0$, we estimated the amount based on literature that indicated the maximum expression is 7-10 times than the basal expression, so we divided maximal expression (≈ 16000μM) by 8 to get the basal expression (2000 μM).

For $\text K_{\text{d1}}$, the equilibrium constant for malic acid and mleR, we used a molecular docking approach to effectively estimate the parameter. By using in silico protein-ligand docking tools, we could determine the binding affinity of two molecules. According to the equation for calculating Gibbs free energy, ΔG° = -RTlnK, we would be able to calculate the equilibrium constant of Malic acid - mleR reaction under a designated temperature.

We chose Pyrx as our in-silico docking tool, due to its easy access and usage (SourceForge). We utilized Benchling’s integrated AlphaFold 2, an in silico protein folding predictor, in order to prepare a PDBQT file of mleR for the docking simulation. AlphaFold 2, from Google’s Deep Mind, is a highly accurate protein prediction software using only the DNA sequence of the protein (Skolnick et al, 2021). We utilized it to Then, the sequence of mleR is imported to AlphaFold in order to produce a hypothetical prediction of the protein’s structure. The PDB file of malic acid is downloaded from the Zinc20 database. Both the files are then imported into Pyrx, with malic acid assigned as a ligand and mleR as a macromolecule.

The docking grid box (a tool for assigning docking area) is set to cover the whole mleR macromolecule. Conventionally, the grid box would only cover specific active sites of the macromolecule that the ligand is designated to bind on, however, due to the lack of literature on the protein complex, we were unable to designate the active site on mleR for malic acid to bind to. For that reason, we allowed malic acid to freely bind to mleR, and used the data for the docking site in which the binding affinity is the lowest (the lower the binding affinity, the more likely the complex is going to be structured that way).
The minimum binding affinity is calculated to be -4.4 kcal/mol, and the equilibrium constant is 0.9926 when the temperature is 26 degrees Celsius.
We assembled all of our parameters into a list:

Using the parameters and the equation, we created a graph comparing the concentration of CcdA and CcdB over the concentration of malic acid.

Conclusion/Discussion

From the model, we can see that at around 1.079 μM of malic acid, the concentration of CcdA overtakes the concentration of CcdB. This indicates our root-localization kill-switch will work as long as the malic acid concentration is above 1.079 μM, as the concentration of CcdA only has to equal the concentration of CcdB to negate any toxicity (Afif et al, 2001).

To estimate the amount of malic acid in the soil, we compared the extent of chemotaxis induced by the banana root exudate and different concentrations of malic acid and we were able to approximate the concentration of malic acid in banana root exudate at 0.279e^0.0627•70 μM, roughly equal to 22.5 μM. Looking back at the graph, we see that our bacteria will survive at this concentration of malic acid, with a very wide margin. The margin provides a wide range of conditions for our bacteria to survive and thrive in without needing to be reapplied with every watering. The malic acid concentration is also near the peak of the curve, producing more CcdA than necessary.

If our bacteria were placed in an area with less than 1.079 μM of malic acid, then the CcdB will slowly kill off the colony and prevent it from spreading further. Our model shows that a root-localization kill-switch based on malic acid is viable as it keeps the bacteria alive near the roots, and kills them in areas without malic acid.

In addition, since the design is based on the malic acid concentration, we were advised to be careful with implementation, as malic acid is a commonly found acid in plants and soil, which may allow our bacteria to survive where it shouldn’t. Our bacteria’s malic acid minimump0u is 1.079 μM of malic acid; compared to other plants, we see that it does not survive with chickpea (0.002 μM), tomatoes (0.14 μM), and wheat (0.25 μM) to name a few. Under our estimations, the basal CcdA production is high enough to kill the bacteria in plants.

With these predictions, we came up with experimental trials for the confirmation of the parameters and viability of the bacteria, but unfortunately, we were unable to finish conducting all trials required for data fitting. Using data fitting, we would have been able to compare the experimental data against literature/estimated data.

Reference

Adeleke, R., Nwangburuka, C., & Oboirien, B. (2017). Origins, roles and fate of organic acids in soils: A review. South African Journal of Botany, 108, 393–406. Link to Source

Dallakian, S. (n.d.). Publications Using PyRx. SourceForge. Link to Source

Gelens, L., Hill, L., Vandervelde, A., Danckaert, J., & Loris, R. (2013). A General Model for Toxin-Antitoxin Module Dynamics Can Explain Persister Cell Formation in E. coli. PLOS Computational Biology, 9(8), e1003190. Link to Source

Renault, P., Gaillardin, C., & Heslot, H. (1989). Product of the Lactococcus lactis gene required for malolactic fermentation is homologous to a family of positive regulators. Journal of Bacteriology, 171(6), 3108–3114. Link to Source

Skolnick, J., Gao, M., Zhou, H., & Singh, S. (2021). AlphaFold 2: Why It Works and Its Implications for Understanding the Relationships of Protein Sequence, Structure, and Function. Journal of Chemical Information and Modeling, 61(10), 4827–4831. Link to Source