Team:Waterloo/Model
Mathematics and Modeling
Environmental Model
An important goal of our project is to determine whether the Linuron concentration in the vicinity of the root nodule is significant enough to damage the Rhizosphere and inhibit nodule growth as plants continue to develop.
In order to model this, we introduce a diffusion-reaction-advection model, which is standard in environmental engineering. A byproduct of this model is that we are able to develop heuristics for the magnitude of Linuron and 3,4-DCA runoff. The model is derived from a similar model by Owsianiak et al [1] for Linuron degradation in bioaugmentation beads. Since this paper also modelled Variovorax sp. [4], we were able to obtain many of our parameters from this work of Owsiniak et al.
Introducing the Model
For completeness, we introduce this model from the ground up. First, we model diffusion through the soil with Fick's Laws of diffusion.
$$ C = [\textrm{Linuron}] $$
$$ D = \textrm{diag}(D_x,D_y,D_z) $$
$$ \nabla\cdot(D\nabla C) = \frac{\partial C}{\partial t} $$
The matrix D is called the diffusion matrix, which determines the speed of diffusion in different directions. Next, we add an advection term to model the transport of Linuron by water.
$$ \nabla \cdot (D\nabla C) - \vec{v}_e \cdot \nabla C = \frac{\partial C}{\partial t}$$
This advection term relies on the hydraulic head, which in our case is a measure of water pressure in the soil. However, water tends to flow around the root nodules, so Linuron transport is diffusion-dominated inside the nodule. Therefore we must use the Navier-Stokes equation in order to find the advection velocity. In order to model degradation, we couple this with our kinetic model to derive this coupled system of PDEs:
$$ \nabla \cdot (D\nabla C) - \vec{v}_e \cdot \nabla C = \frac{\partial C}{\partial t} + \chi_n \frac{V_{max}C}{K_m-C} $$
$$ \varrho \left( \frac{\partial v_e}{\partial t} - v_e \cdot \nabla v_e \right) = \nabla \cdot \sigma(v_e,p)+f $$ $$ \nabla \cdot v_e = 0 $$
Assumptions
In the above model, we assumed an averaged rate of water intake - by necessity due to the complexity of the more realistic rain models. This assumption should produce a reasonable approximation since our timescale is on the order of weeks to months, whereas the rate of rainfall is on the order of days.
We also needed to verify that the rate of water transport into the nodule tissue is not limited by physiological factors such as membrane transport. In order to do this, we used the Octanol-Water partition coefficient in order to predict the rate of membrane transport using Overton's Rule [2]:
$$ P=\frac{K_{ow}D}{\ell} $$
By computing the rate of membrane transport, we concluded that the transport of Linuron into plant cells is not membrane-transport limited. In addition, we assume that the background degradation of Linuron by the same pathway is negligible compared to degradation by Variovorax sp. as found by Owsiniak et al..
Numerical Methods
In order to solve this model we first convert the system into its variational form. The Navier-Stokes terms are solved by a midpoint discretization method called Chorin's Method, while the remainder is solved using the backwards Euler method. The following equations are obtained by integrating the above PDEs and applying integration by parts:
$$\int_\Omega \frac{\partial C}{\partial t} w d\tau = \int_\Omega \left[\nabla \cdot (D \nabla C)w - w v_e\cdot \nabla C - \frac{V_{max}C}{K_m-C}\right] d \tau =$$
$$\int_\Omega \left[(-\nabla C)^T D^T \nabla w -wv_e\cdot \nabla C - \frac{V_{max}C}{K_m-C}\right] d \tau + \int_{\partial \Omega} w( \nabla C)\cdot d \vec{s}$$
We now solve this system in Python using a library known as FeNiCS [3], which is an free, open-source library for solving PDEs with Finite Element Methods. Teams interested in using PDE solvers for their iGEM modelling should consider using FeNiCS as a great alternative to expensive proprietary software. Using ParaView, we export animations of our data from FeNiCS. We also used the Python library meshio to convert meshes generated by Gmsh into a format compatible with FeNiCS. As a note to future iGEM teams considering PDE transport models, your discretization must meet the Courant-Freidrichs-Lewy (CFL) condition in order to be stable.
Results
We now analyze a system of two root nodules, spaced 1 centimeter apart. In the soil we start with an initial supply of Linuron at 1e-10mol/cm^3 isolated just above the top nodule. The nodules are each 1 centimeter in radius and are assumed to be saturated with B. japonicum. Using the Puh-B Parameters specified in the next section, we arrive at the following result for the lower nodule
The following are videos of each of the concentrations over time:
Conclusions
What we can draw from these results is that at the optimized NAT1 concentration the system is very stable, and that the 3,4-DCA concentration falls well below toxicity (see the MIC page). This tells us that implementing the system would not result in a significant pile-up of 3,4-DCA in the soil, which is important information for any future projects since it would aid the optimization of the system.LibA-NAT1 Pathway
We were initially considering three enzymes for Linuron degradation into 3,4-DCA: LibA, PuhA, PuhB, and 2 enzymes for 3,4-DCA degradation: NAT1, NAT2. From preliminary kinetic models, we determined that the use of NAT2 would result in an accumulation of 3,4-DCA that would almost certainly be toxic, and that a PuhA-NAT1 system would not degrade Linuron into 3,4-DCAA at a rate fast enough to mitigate the toxicity of Linuron. The dynamics of the PuhB-NAT1 and LibA-NAT1 systems under fs our preliminary models were identical, and the LibA-NAT1 system was chosen. LibA facilitates the degradation of Liunuron into 3,4-DCA through the following reaction:
It is known that purified LibA is a monomeric linuron hydrolase of ∼55 kDa with a Km and a Vmax for linuron of 5.8 μM and 0.16 nmol per minute for
The second step in the reaction involving NAT1 is characterized in Rodrigues-Lima [5]. The rate for 3,4-dichloroaniline degradation with NAT1 is determined to be 68 ± 8 nmol. per min. per mg. of enzyme.
We assume that the LibA concentration and the NAT1 concentration are equal since they are placed under the same promoter in our system. This choice was done strictly due to time constraints, and our models did suggest that placing LibA and NAT1 under different promoters would result in a more effective system. In order to mitigate the problems of using just a Km value for LibA (since none of the Michaelis-Menton parameters have been found), which would be a blatant disregard for enzyme concentration, we decided to use the Kcat from PuhB since its structure and function are very similar to NAT1.
Computations
We use Julia and the DiffEqBiological.jl package in order to determine the dynamics of our reaction network. Under standard assumptions in biochemical modelling, we reduce our model to the following set of equations for our system with no influx of Linuron:
$$\frac{d[\text{Linuron}]}{dt} = - \frac{V_\max^{\text{Lin}} \cdot [\text{Linuron}]}{K_m^{\text{Lin}} + [\text{Linuron}]} [\text{Linuron}],$$ $$\frac{d[\text{3,4-DCA}]}{dt} = - \frac{V_\max^{\text{DCA}} \cdot [\text{3,4-DCA}]}{K_m^{\text{DCA}} + [\text{3,4-DCA}]} [\text{3,4-DCA}] - k_{\text{DCA}}[\text{3,4-DCA}],$$ $$\frac{d[\text{DCAA}]}{dt} = k_{\text{DCAA}}[\text{DCAA}],$$ where DCAA is dichloroacetanilide, which is what NAT1 degrades 3,4-DCA into. In the case of with influx, the only difference is an influx term in the first equation. With 300 pmol of LibA and NAT1, as well as 600 nmol of Linuron, we see that
where X is the Linuron concentration, Y is the 3,4-DCA concentration, and Z is the 3,4-DCAA concentration.
Under a high influx of Linuron, we see that
and with influx:
Conclusions
Under reasonable amounts of uncertainty in our model, the dynamics of our system were not drastically different. From our model, we determined that in order to engineer a system for the real world, we would need to place LibA and NAT1 under different promoters, ideally with NAT1 under a relatively strong promoter when compared to LibA. In addition, we can conclude that the LibA-NAT1 system could be engineered to degrade Linuron into 3,4-dichloroacetanilide under a wide range of thresholds for toxic concentrations of Linuron and 3,4-DCA. Under the above assumptions of our system, we have demonstrated that a LibA-NAT1 system can efficiently and effectively reduce Linuron concentration in the soil while preserving the rhizosphere. Although we know the ratio of NAT1 and LibA that is required to maintain a strong culture of rhizobia, we were unable to engineer it due to the complexity of the construction.
References
- [1] Owsianiak, M., Dechesne, A., Binning, P. J., Chambon, J. C., Sørensen, S. R., & Smets, B. F. (2010). Evaluation of Bioaugmentation with Entrapped Degrading Cells as a Soil Remediation Technology. Environmental Science & Technology, 44(19), 7622–7627. doi: 10.1021/es101160u
- [2] Grime, J. M. A., Edwards, M. A., Rudd, N. C., & Unwin, P. R. (2008). Quantitative visualization of passive transport across bilayer lipid membranes. Proceedings of the National Academy of Sciences, 105(38), 14277–14282. doi: 10.1073/pnas.0803720105
- [3] Zakharov, P. E. (2018). The FEniCS project. Spark. doi: 10.1515/spark.18.13
- [4] Bers, K., Leroy, B., Breugelmans, P., Albers, P., Lavigne, R., Sørensen, S. R., ... Springael, D. (2011). A Novel Hydrolase Identified by Genomic-Proteomic Analysis of Phenylurea Herbicide Mineralization by Variovorax sp. Strain SRS16. Applied and Environmental Microbiology, 77(24), 8754–8764. doi: 10.1128/aem.06162-11
- [5] Rodrigues-Lima, Fernando & Dairou, Julien & Diaz, Clara & Rubio, Maria & Sim, Edith & Spaink, Herman & Dupret, Jean-Marie. (2006). Cloning, functional expression and characterization of Mesorhizobium loti arylamine N-acetyltransferases: Rhizobial symbiosis supplies leguminous plants with the xenobiotic N-acetylation pathway. Molecular microbiology. 60. 505-12. 10.1111/j.1365-2958.2006.05114.x.