logo

Our mission is to ensure the generation of accurate and precise findings.

Contact Us No. 111, Ren'ai Road, Dushu Lake Science and Education Innovation Zone, Suzhou Industrial Park iGEM@xjtlu.edu.cn
Model Group Oct 3, 2022

Model

Overview

The Model session falls into four diverse but complementary categories, which assist wet lab work and overall insightful optimization of our biofilm system for industrial water purification and heavy metal absorption. We have achieved full-system simulations ranging from molecules to pathways, to populational bacterium working limit, and finally to the mathematical hydrodynamic simulation of the utilized device.

Part 1: Molecular Model

Abstract

Metal binding peptide and Metallothioneins are functional and novel parts of our biofilm system design, in which we aim to comprehensively build and test. An atomic level simulation is tremendously helpful and typically generates substantial insight about how the target component functions, simultaneously proving the feasibility and robustness of the system.

Method & Results
Molecular Docking

Here we present five groups of docking simulation under the computational framework-Autodock vina[1], proving the binding affinity between target heavy metal ions and protein complex respectively.

Receptor Ligand
(1) ~shMT3 Cd2+\stackrel{}{^{}_{}Cd^{2+}_{}}
(2)~shMT3 Zn2+\stackrel{}{^{}_{}Zn^{2+}_{}}
(3)~shMT3 Cu2+\stackrel{}{^{}_{}Cu^{2+}_{}}
(4)~shMT3 Cd2+\stackrel{}{^{}_{}Cd^{2+}_{}} , Zn2+\stackrel{}{^{}_{}Zn^{2+}_{}},Cu2+\stackrel{}{^{}_{}Cu^{2+}_{}}
(5)~shMT3-linker-MBP3 Cd2+\stackrel{}{^{}_{}Cd^{2+}_{}} , Zn2+\stackrel{}{^{}_{}Zn^{2+}_{}},Cu2+\stackrel{}{^{}_{}Cu^{2+}_{}},Ag+\stackrel{}{^{}_{}Ag^{+}_{}}


Number of distinct conformational clusters found = 5, out of 50 runs
Outputting structurally similar clusters, ranked in order of increasing energy.

Figure1.1 shows the binding energy between shMT3 and three heavy metal ions correspondingly.Subsequently, Figure1.2 shows the performance of our designed complex in terms of bioremediation capability.

In conclusion, analysis of the .dlg file from the docking results has successfully proved the reasonable binding affinity of our bioremediation system, while inspiring alternative approach to conduct research on further dynamic docking mechanism.

View bach file/procedure/source code/detailed results

Molecular Dynamics

As part of the metallothionein protein family, shMT3 is Cysteine abundant,from which the existence of Cys-xx-Cys motif support the prediction that sulfhydryl of Cysteine serves as leading factors in terms of chelating Cadmium,Copper,Zinc ions,while two sulfhydryl from different cysteines could be oxidized by oxidants and form a disulfide bond, leading to the defection of heavy ions chelation which could be one possible reason to the lost of function.Alternatively,as shown in the aforementioned docking results, Glutamines are possible binding sites in shMT3 since oxygen atom is more electronegative than sulfur atom. Therefore, it is crucial to comprehensively analyze and optimize its bioremediation performance and relevant insights using computational framework for Molecular Dynamics simulation-Gromacs[2].

General workflow

The increase in the RMSD plot with time shows the protein steadily deviates from its original conformation.The one prominent peak with two ambiguous ones in the histogram suggests the presence of three main conformations which are accessed during the trajectory,possibly indicating the functional motif for chelation with diverse specificity.Higher RMSF values are loop regions with more conformational flexibility, where the structure is not as well defined.This allows a link with experimental spectroscopic techniques which detect the secondary structure of a protein for further functional validation.

After completing the workflow(setup, solvation, minimization, equilibration(NPT,NVT)), we have successfully produced a trajectory (the xtc file) which describes the atomic motion of the system(View detailed procedure of this simulation and guide).From the RMSD,RMSF of the simulation, we could know that the target protein shMT3 will dynamically relax when it be with water molecules under the condition of our project settings, domains which interact with heavy metal ions could gradually be exposed. In this way, the steric hindrance of cadmium ions entering the target can be reduced and the proposed functional motif could assist better chelation performance.

Principal Component Analysis (PCA)

We constructed the covariance matrix of residual movement trajectory of our target functional shMT3 protein using gmx covar and R Bio3D package for visualization.

Cross-correlation analysis was carried out to probe the internal dynamics of our bioremediation system with the results depicted in Figure1.2.3. The positive regions (blue) indicated the strongly correlated motions of residues, while the negative regions (pink) were associated with the anti-correlated movements. The residue movement generally pulls the surrounding atoms in the same direction, so the diagonal regions show highly positively correlated structural motion. For the Metallithonein and Metal Binding peptide we focused on above, the phenomenon was significantly pronounced that free protein displayed the strongest correlated motions, indicating that the correlated motions were reduced when shMT3 chelate with heavy metal ions.


PCA shows the statistically meaningful conformations in the shMT3 trajectory. The principal motions within the trajectory and the vital motions needed for conformational changes can be identified. Two distinct groupings along the PC1 plane, indicating a non-periodic conformational change, are identified. The groupings along the PC2 and PC3 planes do not completely cluster separately, implying that these global motions are periodic. The PC1 is linked to an active site motion that limits the motion to a main chelation interaction.

Gibbs free energy landscape shows the probability density distribution of various conformations of our target functional protein in the direction of PC1 and PC2. As illustrated in the graph, three functional domains are distinct for heavy metal absorption, in consistency with wet lab experiment validation.

View MD bach file/procedure/source code/detailed results

Conclusion & Prospective

In Part1 Molecular level modeling,computational speculation(Autodock,Gromacs,etc) of our designed protein provides theoretical basis and practical guidance for our heavy metal absorption system. We will further proceed to rational design and more intelligent approach like using Deep Learning framework to modify simulation and scoring algorithm based on domain insights are expected, in collaboration with experimental research to construct efficient bioremediation system.

Part 2: Cellular model

View Code

This part focuses on the functions of cells in the CsgA·MBP/MT production and the absorption of metal ions in preparation for the simulation of the absorption ability of the water-purified device.

The production and delivery of CsgA·MBP and CsgA·MT

In our design, only when MBP/MT is out of the cell, they can function as the absorbing units. Due to the equipment and time constraints, we need to model the process of the production and delivery of CsgA·MBP/MT to estimate the concentration of these protein which can be used to adsorb the metal ion and also find the optimal value of input concentration of IPTG.


Fig.2.1 Abstract of the production and delivery of the CsgA·MBP/MT

In our construction, we use tac(lac) promoter to control the production of CsgA·MBP and CsgA·MT. Tac promoter is hybridized by trp promoter and lac promoter[3]. The operon of the tac promoter is almost the same with that of lac promoter. Therefore, the tac promoter can be inhibited by the binding of LacI in the Escherichia coli. In order to adjust the transcription activity of tac promoter, IPTG(Isopropyl β-D-Thiogalactoside) is used to bind with LacI so that LacI will not limit the function of tac promoter.

Aiming to build the model of this process, we makes some assumptions:
1, The model is used to simulate the production and transport of CsgA·MBP/MT.
2, Since the situation of cells in the device is defaulted to be in the stationary phase, the unstable growth rate of E.coli is not considered in these models.
3, LacI can either binds with n1n_1 molecules of IPTG to form the complex(LacIn1IPTGLacI·n_1 IPTG) or unbound so that intermediate states where less than n1n_1 molecules are bound can be neglected[4].
4, The transporter can export CsgAMBP/MTCsgA·MBP/MT due to the existence of CsgA.
5, The model ignores the detailed step that CsgB nucleates CsgA and assumes that CsgA will be continuously exported by transporters, so the export rate of CsgA·MBP and CsgA·MT is directly proportional to the concentration of CsgA·MBP/MT in the cells.
6,The model neglects the diffusion of IPTG into and out of cells, but put some adjustments on the values of input concentration of IPTG.

Based on the above mechanism and assumptions, two equations are combined to model the process that inducer binds repressor to prevent the inhibition on tac promoter. The first equation is to describe the relationship between IPTG and LacI. According to the third assumption, Hill function can be used to determine active LacI which can bind with tac promoter.


Fig.2.2 The equation of the concentration active LacI

In this equation, [LacI][LacI^*] is the concentration of active LacI,[LacI]t[LacI]_t represents the total concentration of LacI (including LacI binding with IPTG and active LacI), [IPTG][IPTG] is the concentration of IPTG input, nn is the hill coefficient and KR1K_{R1} is the half-saturation of IPTG.

The other equation denotes the transcription by the tac promoter. Although the reaction of the binding between LacI and tac promoter is not the process of enzyme kinetics, Michaelis-Menten equation can be used to describe the repression on tac promoter.


Fig.2.3 The equation of promoter activity

β\beta represents the maximal transcription rate of tac promoter and KX1K_{X1} is the half-saturation of active LacI.

By combining the two equations, we can understand the relationship between IPTG and tac promoter activity.


Fig.2.4 The equation of promoter activity combined with the two above equations

Here, we can predict the optimal concentration of input IPTG.


Fig.2.5 The prediction of the promoter activity under different concentrations of IPTG

As shown in the Fig.2.5, with the increase of the concentration of IPTG input, the promoter activity grows up but reaches the peak after [IPTG] = 1 mM (approximately) and remians maximum after 1 mM. Since too high concentration of IPTG will affect the growth of cells which cannot be expressed in this model, the optimal concentration of IPTG needs to cause the promoter activity at maximum and is not too high. Therefore, 1 mM may be a suitable concentration of IPTG input.(Experiment group had experiments related to it in the page).

The following equations are used in the production and secretion of CsgA·MBP/MT

Equations Description
d[mRNAcsgAMBP]dt=β11+[LacI]tKx111+([IPTG]KR1)n1r1[mRNAcsgAMBP]\frac{d[mRNA·CsgA·MBP]}{dt} = \frac{\beta_1}{1+\frac{[LacI]_t}{K_{x1}}·\frac{1}{1+(\frac{[IPTG]}{K_{R1}})^{n_1}}} - r_1[mRNA·CsgA·MBP] Change of the concentration of the mRNA·CsgA·MBP based on IPTG binding LacI to reduce the binding between LacI and tac promoter according to Hill Equation and Michaelis–Menten Equation. β1\beta_1 represents the maximal transcription rate of tac promoter.[LacI]t[LacI]_t is the normal concentration of LacI in one Escherichia coli. Kx1K_{x1} is the unbinding LacI concentration at which transcription rate of tac promoter is at half-maximum. n1n_1 is the hill coefficient of IPTG and LacI. KR1K_{R1} stands for IPTG concentration at which the concentration of active LacI is at half-maximum. r1r_1 is the degradation rate of mRNAcsgAMBPmRNA·CsgA·MBP
d[csgAMBP]dt=k1[mRNAcsgAMBP]r2[csgAMBP]\frac{d[CsgA·MBP]}{dt} = k_1[mRNA·CsgA·MBP] - r_2[CsgA·MBP] Change of concentration of csgAMBPCsgA·MBP due to translation and degradation. k1k_1 is the translation rate of csgAMBPCsgA·MBP. r2r_2 stands for the degradation rate of csgAMBPCsgA·MBP
d[mRNAcsgAMT]dt=β11+[LacI]tKx111+([IPTG]KR1)n1r3[mRNAcsgAMT]\frac{d[mRNA·CsgA·MT]}{dt} = \frac{\beta_1}{1+\frac{[LacI]_t}{K_{x1}}·\frac{1}{1+(\frac{[IPTG]}{K_{R1}})^{n_1}}} - r_3[mRNA·CsgA·MT] Change of the concentration of the mRNA·CsgA·MT based on IPTG binding LacI to reduce the binding of LacI and tac promoter according to Hill Equation and action of mass. The process and values of most parameters are similar to that of csgAMBPCsgA·MBP. r3r_3 is the degradation rate of csgAMTCsgA·MT which is the only value different from that of csgAMBPCsgA·MBP.
d[csgAMT]dt=k2[mRNAcsgAMT]r4[csgAMT]\frac{d[CsgA·MT]}{dt} = k_2[mRNA·CsgA·MT] - r_4[CsgA·MT] Change of concentration of csgAMTCsgA·MT due to translation and degradation. k2k_2 is the translation rate of csgAMTCsgA·MT. r4r_4 stands for the degradation rate of csgAMTCsgA·MT
d[csgAMBP/MT]exdt=k3[csgAMBP/MT]\frac{d[CsgA·MBP/MT]_{ex}}{dt} = k_3[CsgA·MBP/MT] Transport of CsgA·MBP/MT into outside of the cell. k3k_3 is the export rate of csgAMBP/MTCsgA·MBP/MT


To be noticed, units of most of the above equations is μMs1\mu M·s^{-1} while the unit of the export rate of CsgA·MBP/MT is μmols1\mu mol·s^{-1}. This is because we standardize the export rate k4k_4. According to law of mass action, the speed of one reaction is connected with the concentration of each substrate because the the higher concentration of substrate is, the higher probability molecules of substrates have to collide and react with each other. This means that the determination of value of the concentration of CsgA·MBP/MT is based on the volume of the container where cells live in actual situations. Therefore, the export rate is standardized to adapt different cases.

These are the parameters we used:

Parameters Description Source
β1=0.0014(μMs1)\beta_1 = 0.0014 (\mu M·s^{-1}) Maximal transcription rate of tac promoter Estimated from [5]
[LacI]t=0.01(μM)[LacI]_t = 0.01 (\mu M) The normal concentration of LacI in one Escherichia coli Estimated from [6]
Kx1=5E5(μM)K_{x1} = 5E-5(\mu M) The unbinding LacI concentration at which transcription rate of tac promoter is at half-maximum Estimated from [5:1]
KR1=1(μM)K_{R1} = 1(\mu M) Dissociation constant of IPTG and LacI Estimated from [7]
n1=2n_1 = 2 Hill coefficient of IPTG and LacI Estimated from [7:1]
r1=4.4E3(s1)r_1 = 4.4E-3(s^{-1}) The degradation rate of mRNAcsgAMBPmRNA·CsgA·MBP Estimated from [8]
k1=0.57(s1)k_1 = 0.57(s^{-1}) The translation rate of mRNAcsgAMBPmRNA·CsgA·MBP Estimated from [8:1]
k2=0.44(s1)k_2 = 0.44(s^{-1}) The translation rate of mRNAcsgAMTmRNA·CsgA·MT Estimated from [8:2]
r2=6.3E5(s1)r_2 = 6.3E-5(s^{-1}) The degradation rate of csgAMBPCsgA·MBP Estimated from [9]
r3=3.4E3(s1)r_3 = 3.4E-3(s^{-1}) The degradation rate of mRNAcsgAMTmRNA·CsgA·MT Estimated from [8:3]
r4=4.88E5(s1)r_4 = 4.88E-5(s^{-1}) The degradation rate of csgAMTCsgA·MT Estimated from [9:1]
k3=1.8E19(Ls1)k_3 = 1.8E-19(L·s^{-1}) The export rate of CsgA·MBP/MT Estimated from Team:TU_Delft (2015)


According to values of parameters and equations in this model, we can predict the amount of substance of CsgA·MBP and CsgA·MT:


Fig.2.6 The changes of amount of substance of CsgA·MBP/MT with time under 1 mM IPTG

According to the Fig.2.6, with certain time period and the size of contianer, the model can estimate the concentration of CsgA·MBP and CsgA·MT in preparation of fitting the experimetal data in absorption.

The absorption of silver ion

In this part, we build a model about the absorption of silver ions and verify this model according to the fitting of the experimental data. The values of key parameters acquired from fitting lay a foundation for further development in simulation of water-purified device.

The mechanism of the absorption implemented by MBP is that firstly MBP changes silver ion into AgNP (This process is irreversible) and then absorbs the AgNP(This process is reversible)[10].


Fig.2.7 The biological pathway of the absorption of AgNP

Since the principle of the reaction mainly results from that the tyrosine in the MBP transfer its electron to the silver ion, we use [CsgAMBP][CsgA·MBP] to represent the concentration of CsgAMBPCsgA·MBP which does not react with silver ions and [CsgAMBP+][CsgA·MBP^+] stands for the concentration of CsgAMBPCsgA·MBP which has been reacted with silver ions. Moreover, the part in the CsgAMBPCsgA·MBP functioning as adsorption is different from the part functioning as reaction in the protein, so we use [CM][CM] to stand for the concentration of CsgA·MBP which does not absorb silver nanoparticles.

However, in the process of building of the models, we encounter some troubles. From experiment, we can only know how many silver ions are absorbed by MBP (whether all these silver ion become AgNP can not be known). The other difficulty is that we are unable to acquire the values of key parameters about the reaction and absorption from previous articles. As a result, in our model, we combine the two mechanisms into one absorption process.


Fig.2.8 The modified biological pathway of the absorption of AgNP

Here, [CsgAMBP][CsgA·MBP] stands for CsgA·MBP which does not chelate silver ions or AgNP. [CsgAMBPn2AgNP][CsgA·MBP·n_2 AgNP] represents the CsgA·MBP which chelates silver ions or AgNP.

Considering the simplification and the whole process, assumptions are made:
1, For the process of absorption, we neglect the reaction which changes silver ions into silver nanoparticles and mix it into the reaction of absorption.
2, CsgAMBPCsgA·MBP can either binds with n2n_2 molecules of AgNP(silver ions) to form the complex(CsgAMBPn2AgNP/Ag+CsgA·MBP·n_2 AgNP/Ag^+) or unbound so that intermediate states where less than n2n_2 molecules are bound can be neglected[4:1].
3, When we consider the irreversible absorption, we do not matter whether the units absorbed are AgNPs or silver ions. The binding affinity is used to describe both of the objectives.
4,The concentration of AgNP in our model is more likely to be referred to the concentration of Ag(0)(simple substrate of silver) rather than how many nanoparticles(whose dimeter is larger than 1 nm) produced in the system, as it is difficult to obtain the number of nanoparticles.
5,This model do not consider the process that cells themselves transport silver ions into their bodies.

The following equations are used in the absorption of silver ions

Equations Description
d[CsgAMBPn2AgNP]dt=kon1[CsgAMBP][Ag+]n2koff1[CsgAMBPn2AgNP]\frac{d[CsgA·MBP·n_2 AgNP]}{dt} = k_{on1}[CsgA·MBP][Ag^{+}]^{n_2} - k_{off1}[CsgA·MBP·n_2 AgNP] The production of mixture of CsgA·MBP and AgNP(Ag+Ag^{+}). kon1k_{on1} represents the reaction rate of association between CsgA·MBP and AgNp(Ag+Ag^{+}). koff1k_{off1} stands for the reaction rate of dissociation between CsgA·MBP and AgNp(Ag+Ag^{+}). n2n_2 is the hill coefficient of CsgA·MBP and AgNP(Ag+Ag^{+})
d[AgNPabsorbed]dt=n2(kon1[CsgAMBP][Ag+]n2koff1[CsgAMBPn2AgNP])\frac{d[AgNP·absorbed]}{dt} = n_2(k_{on1}[CsgA·MBP][Ag^{+}]^{n_2} - k_{off1}[CsgA·MBP·n_2 AgNP]) The change of the concentration of the AgNP(Ag+Ag^{+}) absorbed. The values of parameters are the same with the first equation
d[CsgAMBP]dt=(kon1[CsgAMBP][Ag+]n2koff1[CsgAMBPn2AgNP])\frac{d[CsgA·MBP]}{dt} = -(k_{on1}[CsgA·MBP][Ag^{+}]^{n_2} - k_{off1}[CsgA·MBP·n_2 AgNP]) The change of the concentration of CsgA·MBP which does not binds with AgNP(Ag+Ag^{+}). The values of parameters are the same with the first equation


Since it is impractical for us to measure precisely the changes of the concentrations of silver ions and CsgA·MBP/MT(In other words, we cannot acquire the values of kon1k_{on1} and koff1k_{off1} directly). Instead, we can use binding affinity to alternative the process.

In this part, we replace kon1\pmb{k_{on1}} and koff1\pmb{k_{off1}} with KD1\pmb{K_{D1}} which is calculated by koff1\pmb{k_{off1}} dividing kon1\pmb{k_{on1}}[11]. According to the equation: d[AgNPabsorbed]dt=n2(kon1[CsgAMBP][Ag+]n2koff1[CsgAMBPn2AgNP])\pmb{\frac{d[AgNP·absorbed]}{dt} = n_2(k_{on1}[CsgA·MBP][Ag^{+}]^{n_2} - k_{off1}[CsgA·MBP·n_2 AgNP])}, when d[AgNPabsorbed]dt\frac{d[AgNP·absorbed]}{dt} becomes zero (which means the whole system reaches steady state and cannot absorb AgNp(Ag+\pmb{Ag^+}) any more), we can obtain that [AgNPabsorbed]=n2[CsgAMBP][Ag+]n2KD1+[Ag+]n2\pmb{[AgNP·absorbed] = \frac{n_2[CsgA·MBP][Ag^+]^{n_2}}{K_{D1} + [Ag^+]^{n_2}}} which represents the maximal concentration of silver ions absorption. Here, [Ag+][Ag^+] represents the concentration of silver ions which is unbound. This equation is used to fit with experiment data to obtain the value of the binding affinity and hill coefficient.


Fig.2.9 The equation for the maximal absorption of silver ions

According to the results from experiment, we use matlab to fit the data by this model:

Fig.2.10 The fitting of experimental data

According to the fitting, we can obtian the significant values of binding affinity and hill coefficient in the absorption of silver ions.

These are the parameters we used and values acquired from fitting:

Parameters Description Source
KD1=6.655(μM)K_{D1} = 6.655(\mu M) Binding affinity between CsgA·MBP and AgNP(Ag+Ag^+) Estimated from experiment
n2=0.5086n_2 = 0.5086 Hill coefficient between CsgA·MBP and AgNP(silver ions) Estimated from experiment
The absorption of Cu2+Cu^{2+} and Cd2+Cd^{2+} by CsgA·MT

In this part, we build a model about the absorption of copper and cadmium ions and verify this model according to the fitting of the experimental data. The values of key parameters acquired from fitting lay a foundation for further development in simulation of water-purified device.

The mechanism of the absorption by MT is that it chelates the metal ions(Cu2+,Cd2+Cu^{2+},Cd^{2+}) and the process is reversible[12].


Fig.2.9 The biological pathway for the absorption of Cu2+Cu^{2+}


Fig.2.10 The biological pathway for the absorption of Cd2+Cd^{2+}

The assumptions we made:
1, csgAMTCsgA·MT can either binds with n3n_3/n4n_4 molecules of Cu2+/Cd2+Cu^{2+}/Cd^{2+} to form the complex(csgAMTn3Cu2+/n4Cd2+CsgA·MT·n_3 Cu^{2+}/n_4 Cd^{2+}) or unbound so that intermediate states where less than n3/n4n_3/n_4 molecules are bound can be neglected[4:2].
2, This model do not consider the process that cells themselves transport copper and cadmium ions into their bodies.

The following equations is used to describe the chelation process.

Equations Description
d[csgAMTn3Cu2+]dt=kon2[csgAMT][Cu2+]n3koff2[csgAMTn3Cu2+]\frac{d[CsgA·MT·n_3 Cu^{2+}]}{dt} = k_{on2}[CsgA·MT][Cu^{2+}]^{n_3} - k_{off2}[CsgA·MT·n_3 Cu{2+}] The production of mixture of CsgA·MT and Cu^{2+}. kon2k_{on2} represents the reaction rate of association for CsgA·MT and Cu^{2+}. koff2k_{off2} stands for the reaction rate of dissociation for CsgA·MT and Cu^{2+}. n3n_3 is the hill coefficient of CsgA·MT and Cu^{2+}
d[Cu2+absorbed]dt=n3(kon2[csgAMT][Cu2+]n3koff2[csgAMTn3Cu2+])\frac{d[Cu^{2+}·absorbed]}{dt} = n_3(k_{on2}[CsgA·MT][Cu^{2+}]^{n_3} - k_{off2}[CsgA·MT·n_3 Cu^{2+}]) The change of the concentration of the Cu^{2+} absorbed. The values of parameters are the same with the first equation for the chelation of Cu^{2+}
d[csgAMT]dt=(kon2[csgAMT][Cu2+]n3koff2[csgAMTn3Cu2+])\frac{d[CsgA·MT]}{dt} = -(k_{on2}[CsgA·MT][Cu^{2+}]^{n_3} - k_{off2}[CsgA·MT·n_3 Cu^{2+}]) The change of the concentration of CsgA·MT which does not chelate with Cu^{2+}. The values of parameters are the same with the first equation for the chelation of Cu^{2+}
d[csgAMTn4Cd2+]dt=kon3[csgAMT][Cd2+]n4koff3[csgAMTn4Cd2+]\frac{d[CsgA·MT·n_4 Cd^{2+}]}{dt} = k_{on3}[CsgA·MT][Cd^{2+}]^{n_4} - k_{off3}[CsgA·MT·n_4 Cd{2+}] The production of mixture of CsgA·MT and Cd^{2+}. kon3k_{on3} represents the reaction rate of association for CsgA·MT and Cd^{2+}. koff3k_{off3} stands for the reaction rate of dissociation for CsgA·MT and Cd^{2+}. n4n_4 is the hill coefficient of CsgA·MT and Cd^{2+}
d[Cd2+absorbed]dt=n4(kon3[csgAMT][Cd2+]n4koff3[csgAMTn4Cd2+])\frac{d[Cd^{2+}·absorbed]}{dt} = n_4(k_{on3}[CsgA·MT][Cd^{2+}]^{n_4} - k_{off3}[CsgA·MT·n_4 Cd^{2+}]) The change of the concentration of the Cd^{2+} absorbed. The values of parameters are the same with the first equation for the chelation of Cd^{2+}
d[csgAMT]dt=(kon3[csgAMT][Cd2+]n4koff3[csgAMTn4Cd2+])\frac{d[CsgA·MT]}{dt} = -(k_{on3}[CsgA·MT][Cd^{2+}]^{n_4} - k_{off3}[CsgA·MT·n_4 Cd^{2+}]) The change of the concentration of CsgA·MT which does not chelate with Cd^{2+}. The values of parameters are the same with the first equation for the chelation of Cd^{2+}


Just like how we solve with the data about the absorption of AgNP(Ag+Ag^{+}), we can recognize that the maximal concentration of absorbed Cu2+\pmb{Cu^{2+}} and Cd2+\pmb{Cd^{2+}} is [Cu2+absorbed]=n3[CsgAMT][Cu2+]n3KD2+[Cu2+]n3\pmb{[Cu^{2+}·absorbed] = \frac{n_3[CsgA·MT][Cu^{2+}]^{n_3}}{K_{D2} + [Cu^{2+}]^{n_3}}} and [Cd2+absorbed]=n4[CsgAMT][Cd2+]n4KD3+[Cd2+]n4\pmb{[Cd^{2+}·absorbed] = \frac{n_4[CsgA·MT][Cd^{2+}]^{n_4}}{K_{D3} + [Cd^{2+}]^{n_4}}}.[Cu2+][Cu^{2+}] and [Cd2+][Cd^{2+}] are concentrations of unbound copper ions and cadmium ions. Combined with experiment data, we can get the values of KD1\pmb{K_{D1}}, KD2\pmb{K_{D2}}, n3\pmb{n_3} and n4\pmb{n_4}.

Fig.2.11 The maximal absorption for copper ions

Fig.2.12 The maximal absorption for cadmium ions

These are the parameters we used:

Parameters Description Source
KD2=0.3559(μM)K_{D2} = 0.3559(\mu M) Binding affinity between CsgA·MT and Cu2+Cu^{2+} Estimated from experiment
KD3=0.1676(μM)K_{D3} = 0.1676(\mu M) Binding affinity between CsgA·MT and Cd2+Cd^{2+} Estimated from experiment
n3=0.4614n_3 = 0.4614 Hill coefficient in the binding of between Cu2+Cu^{2+} and CsgA·MT Estimated from experiment
n4=0.3182n_4 = 0.3182 Hill coefficient in the binding of between Cd2+Cd^{2+} and CsgA·MT Estimated from experiment

Part 3 : Population Survival Model
>

View Code

Introduction

In industrial wastewater environments, metal ions have inhibitory effects on the growth of E. coli. However, the degree of tolerance to metal ions in our genetically modified E. coli still needs to be explored. Here, we first developed a model to fit the experimental data to understand the growth situation of E. coli parametrically. Then, we performed a dependence analysis of the parameters in the model and considered how these parameters change with increasing metal ions concentrations. In addition, a range of E. coli growth curves was predicted. The biological processes were parametrically evaluated in this part, leading to a more in-depth interpretation of the experimental results and further predictions to guide implementation decisions in the real world.

Model Development

Consider the functional relationship between cell number (denoted by n(t)) and time (denoted by t). In the case of unlimited resources, the growth rate r is constant, resulting in an exponential growth of E. coli. However, in wastewater environments, r depends on the cell density and varies with time. r=r(n)=r0r = r(n) = r_0×(1- nK\frac{n}{K}), where r0r_0 is inherent growth rate.

The equation of E. coli growth hence become dndt\frac{dn}{dt} = n(t)×r0n(t)×r_0×(1-n(t)K)\frac{n(t)}{K}). Therefore, the logistic growth curve of E. coli can be shown by n(t)=K(1+(Kn0)n0×er0t)n(t) = \frac{K}{(1+\frac{(K-n_0)}{n_0} ×e^{-r_0 t} )}

Since the cell number is equivalent to the cell density or approximately the OD600OD_{600} readings, the OD600OD_{600} values and the cell density follow the same n(t)n(t) expression [13] .

To assess the extent to which the growth of E. coli was inhibited, we established the SAD model to describe the state change of E. coli after the addition of metal ions[13:1]. It is assumed that all E. coli were knocked into a suppressed state, S, after the addition of metal ions. After that, two possible states of E. coli were available for conversion. Some E. coli regained their activity, woke up and returned to the active state A, and then continued the cycle of multiplying and growing. Some other E. coli were changed from a suppressed state to a dead state D.

Parameters Description
r0r_0 Inherent growth rate
n0n_0 Initial cell density
KK Asymptotic cell density
nsn_s Cell density in the suppressed state
nan_a Cell density in the active state
ndn_d Cell density in the dead state
α\alpha Killing rate, from the suppressed state to the dead state
β\beta Weak-up rate, from the suppressed state to the active state


Then the SAD model can be quantitatively described as follows: [13:2]:

dnsdt=β×nsα×ns\frac{dn_s}{dt} = -β×n_s - α×n_s

dnadt=β×ns+r0×(1na(Knsnd))×na\frac{dn_a}{dt} = β×n_s+r_0×(1-\frac{n_a}{(K-n_s-n_d )}) ×n_a

dnddt=α×ns\frac{dn_d}{dt} = α×n_s

n(t0)=n0n(t_0) = n_0

Fig. 3.1 The three-state transition (SAD) model.

Results and Analysis

Here is the framework for analyzing natural biological phenomena with uncertain outcomes. The real world generates data and is evaluated by our experiments. The field of data processing uses these data to come up with a quantitative description model. Once we develop the model, we use analysis tools to explore deeper relations and evaluate them. Moreover, the results from this analysis lead to predictions and decisions about the real world.

Fig. 3.2 The relation between the real phenomena, data processing and analysis.

Fitting Analysis of Experiment Data

As shown in Fig. 3.3 A, the survival situation of the engineered MG1655 pGEX-4T-1-csgA-MBP3 can be observed:

No significant effect was observed for Ag ions below 5 μM.

At lower concentrations (Ag ions between 10 μM and 20 μM), E. coli was inhibited for the first approximately 3 hours, then growth resumed.

In the first approximately 6 hours (with an initial OD600OD_{600}=0.05), Ag ions at a concentration of 40μM completely inhibit the growth of E. coli.

Under 80 μM Ag ions treatment, E. coli grew again after approximately 13 hours of completely inhibited.

These observations are consistent with most previous works of literature [13:3],[14]. Also, it can be found that after various periods of suppression, the E. coli population eventually began to grow again. Therefore, it is likely that the suppression of the growth of E. coli by Ag ions is only transient under these experimental conditions. This observation also confirms the rationality of the three-state transition (SAD) model.

Furthermore, although the growth and proliferation of bacteria are widely monitored by measuring the optical density (OD) at 600 nm, the E. coli growth situation may still not be genuinely displayed due to the obstacles in the OD method, such as multiple scattering [14:2]. In our experiment, the growth curves of E. coli, especially under high concentrations of Ag ions (above 40 μM), deviated to some extent from the sigmoid shape (Fig. 3.3 A).

To minimize the error and find the most likely value for each time of the experiment, we draw the method mentioned in the literature to overcome the barriers in evaluating bacterial growth with microplate readers[14:2]. The first-order time derivative of the OD values and its correspondence with the model parameters were evaluated. This method is consistent with the traditional sigmoid shape of E. coli growth. Furthermore, it provides a better interpretation of OD growth curves when the experimentally measured growth curves do not follow the accepted sigmoid shape due to multiple scattering. Fig. 3.3 B, C, and D show the growth curves of three kinds of E. coli under different silver ion concentrations fitted by the Gompertz model,A×exp[exp[Um×eK(λt)+1]]+n0A×exp⁡[ -exp⁡[\frac{Um×e}{K} (λ-t)+1] ]+ n_0

Fig. 3.3 Growth curves fitted by Gompertz model. (A) The experiment data of OD600OD_{600} values of MG1655 pGEX-4T-1-csgA-MBP3 at different Ag ions concentration. (B), (C) and (D) Model-fitted growth curves of MG1655 pGEX-4T-1-csgA-MBP3, MG1655 and MG1655 pGEX-4T-1, respectively.

Dependence Analysis of Parameters in Model

To further study the variation of growth curves and evaluate the effect of genetic modification on the tolerance of Ag ions, we selected the lag time (λ) and the maximum growth rate (Um) as parameters to characterize the growth of E. coli. These two parameters are evaluated by fitting analysis in each concentration of these three kinds of E. coli.

The experiments are conducted on three types of E. coli: (a) MG1655: Wide type of E. coli (b) MG1655 pGEX-4T-1: MG1655 imported the empty plasmid (c) MG1655 pGEX-4T-1-csgA-MBP3: Genetically engineered MG1655, which is used to treat wastewater with Ag ions

Focus on MG1655 pGEX-4T-1-csgA-MBP3, as shown in Fig. 3.4, the maximum specific growth rate (Um) is essentially the same for different concentrations of Ag ions in the growth medium. However, as the concentration of Ag ions in the growth medium increased, the lag time (λ) increased significantly [13:4]. A similar variation of these two parameters can be observed in MG1655 and MG1655 pGEX-4T-1.

In addition, compare the parameters’ variation of these three kinds of E. coli. It can be found that the lambda of MG1655 pGEX-4T-1 in each Ag ion concentration is more similar to the lambda value of MG1655. At lower concentrations (Ag ions between 0 μM and 20 μM), the lag time in MG1655 pGEX-4T-1-csgA-MBP3 is more significant. Besides that, in these three kinds of E. coli, the maximum specific growth rate (Um) remained stable (approximately between 1.0-1.2 in experiment conditions).

Therefore, it is reasonable to believe that the operation of import plasmid has a limited effect on the Ag ion tolerance of E. coli. Compared with wild type MG1655, the variation of MG1655 pGEX-4T-1 growth curves was mainly caused by genetic modification. In addition, the inhibition effect caused by Ag-treatment was mainly reflected in the lag time elongation without changing the maximum specific growth rate. Our engineered MG1655 pGEX-4T-1-csgA-MBP3 could adapt and grow in the environment of 80 μM Ag ions after a period of inhibition.

Fig. 3.4 Dependence of the lag time (λ) and the maximum specific growth rate (Um) of E. coli on Ag ions concentration in growth medium

Model Predictions

Based on the three-state transition (Suppressed-Active-Dead) model, it is possible to make several predicted OD600OD_{600} and CFU curves. As explained above, 𝛼 is the killing rate from the suppressed state to the dead state, and 𝛽 is the wake-up rate from suppressed state to the active state. Note: all the E. coli (nsn_s + nan_a + ndn_d) are included in plotting the OD600OD_{600} curves. Only alive E. coli (nsn_s + nan_a) are included in CFU curves.

As we have previously analyzed, the shape of the growth curves is sigmoid. Also, the exponential growth phase of MG1655 pGEX-4T-1-csgA-MBP3 is essentially parallel to the exponential growth phase of MG1655 at all different Ag ions concentrations, namely at different wake-up rates (β). This fact indicates that the maximum specific growth rate is not affected by the inhibition caused by Ag treatment or the wake-up rate (β). Differently, the lag time was longer for smaller wake-up rates (β), which corresponded to higher Ag ion concentrations.

In addition, we also predicted the CFU curves, which can reflect the time-killing situation of Ag ions to E. coli. As shown in Fig. 3.5 B, the CFU counts decreased in the first few hours due to death caused by metal ions and then increased again at roughly the same rate as untreated E. coli, which is consistent with reality. Also, since the killing rate is mainly determined by the types of metal ions, the concentrations of Ag ions, and the size of AgNPs [13:5], it is reasonable to assume a constant killing rate in our experimental conditions. As the wake-up rate (β) decreased, it took longer for E. coli to grow again, meaning it took longer to adapt to the environment (Fig. 5B).

For different concentrations of Ag ions, it is natural to consider the CFU curves at different killing rates (a). As shown in Fig. 3.5 C, with an increasing killing rate, the turning point of the CFU curves was advanced, accompanied by a smaller number of cells at the time of growing again. Fig. 3.5 D magnifies the CFU curves from 0 to 9 hours, which shows the variations more clearly.

Fig. 3.5 Prediction curves. (A) Predicted kinetic growth curves with Ag added in the lag phase (a = 0.77). (B) Predicted CFU counts with various 𝛽 (a = 0.77). (C) and (D) Predicted CFU counts with various a (𝛽 = 0.01).

Conclusion and discussion

To conclude, our model provides a quantified description of the tolerance activity of E. coli to Ag ions. The toxic effect of Ag ions on E. coli fed back mainly on the variation of two parameters, the lag time and the maximum specific growth rate. Model fitting analysis of the experimental data revealed that the lag time of E. coli growth increased significantly with increasing Ag-treatment concentration. In the range of Ag ion concentrations from 0 μM to 80 μM, the E. coli could grow again, and the maximum growth rate remained basically unchanged. Furthermore, The growth situations of MG1655 pGEX-4T-1, which imported the empty plasmid, was much closer to the wild type, indicating that Ag ions mainly caused the growth inhibition of MG1655 pGEX-4T-1-csgA-MBP3, the effect caused by the operation of import plasmid was very limited.

Additionally, it is expected that more parameters characterizing the growth of E. coli can be included in our model, allowing a more detailed exploration of small changes in actual growth and the biological implications reflected. It is also important to emphasize the extensibility of our model for assessing metal ions other than Ag ions. The tolerance activity of E. coli to other heavy metal ions can be similarly obtained. By establishing a network of relationships between multiple heavy metal ion concentrations and the corresponding E. coli survival parameters, we can predict the survival of E. coli in natural wastewater environments and guide practical operations.

Part 4: System Model

View Code

Demonstration of the ability of the absorption

Our goal of our water purified device is to ensure that after implemented, the concentration of metal ions in the water is reduced to what is below that of the national standard. In the laboratory, it is difficult to simulate the situation that continuous water pass through the device. As a result, we design a simple model to simulate this detailed process.

According to previous calculations, the maximal absorption of the metal ions is [metalabsorbed]=n[CsgAMBP/MT][metal]nKD+[metal]n\pmb{[metal·absorbed] = \frac{n [CsgA·MBP/MT][metal]^{n}}{K_D + [metal]^{n}}} where [CsgAMBP/MT][CsgA·MBP/MT] represents the concentration of the CsgA·MBP/MT, [metal][metal] is the concentration of metal ions in the device, nn is the hill coefficient and KDK_D is the binding affinity. Therefore, the rate of change in concentration is [metalabsorbed]=n[CsgAMBP/MT][metal]n(KD+[metal]n)t\pmb{[metal·absorbed] = \frac{n[CsgA·MBP/MT][metal]^{n}}{(K_D + [metal]^{n})*t}}. In this case, t=36/8ht = 36/8 h.

However, the above equation is not enough to simulate the real absorption. When the water with metal ions pass into the device, the concentration will change with the increase of the input position and thus causes the change of the concentration of CsgAMBP/MTCsgA·MBP/MT. As a result, with time passes, the rate of absorption will change with the position and changes in the concentration of metal ions. It is the reason why we need to simulate the distribution of the concentration of metal ions and thus predict the concentration of output metal ions.

In addition, in our design, cells used grow on the MBBR carriers to filter the water which increase the complexity of the absorption situation. In reality, the water-purified equipment is filled with MBBR carriers and these carrier will move with the water when current with metal ions is input. Since the process involves with vortex and the corresponding movement of MBBR carriers, it is difficult to consider the elaborate procedures of the absorption in the device. Therefore, necessary simplifications are required to help with the simulation.

The assumptions we made:
1, we neglect the movement of and the detailed structure of MBBR carriers. Instead, we assume that the device is full of absorption units evenly at first as MBBR carriers occupy almost the whole space.
2, we assume that directions of the currents are the same from the inlet and outlet so that we can neglect the vortex in the device.
3, The diffusion is not considered in the equation as the effect of is limited compared with the speed of flow

Here, we get the equation as below:

Fig.4.1 The equation of the absorption of metal ions in the water-purified device

On the left hand side, the term represents how the concentration of metal ions changes at each point with the pass of time. On the right hand side, the first term describes the change of concentration of metal ions due to flow and the second term is the absorption rate.

With the help of the matlab, we find the change of the concentration of the metal ions in the device and suitable values of the parameters(the values of K\pmb{K} and n\pmb{n} result from experiment).

Fig.4.2 An example of simulation of the distribution of concentration of silver ions in the device

In this example,the end point of the distance(x-axis) represents the output concentration of silver ions in the current which and input concentration of silver ions is random within the range of 0 ~ 12 μM\pmb{\mu M} for each time steps with v = 100 m/h and initial [CsgAMBP/MT][CsgA·MBP/MT] = 1.2M per position. The output concentration proves the absorption ability of our device as it is less than that of national stanard(Here is 0.9 μM\mu M) The values of speed of water and initial concentration of CsgA·MBP/MT can be changed in the equation according to real needs.

Under the guidance of this model, we can determine the suitable input of the number of MBBR carriers, number of cells incubated on these carriers and the current speed to ensure the purification of the output water on the metal ions.

We hope that the method of how we create the model to simulate the output concentration of metal ions can give a good idea for future teams to solve similar problems within absorption in fluid mechanics and water-purified device.

Appendix

View github repository for the Model session

Parameters Description Source
β1=0.0014(μMs1)\beta_1 = 0.0014 (\mu M·s^{-1}) Maximal transcription rate of tac promoter Estimated from [5:2]
[LacI]t=0.01(μM)[LacI]_t = 0.01 (\mu M) The normal concentration of LacI in one Escherichia coli Estimated from [6:1]
Kx1=5E5(μM)K_{x1} = 5E-5(\mu M) The unbinding LacI concentration at which transcription rate of tac promoter is at half-maximum Estimated from [5:3]
KR1=1(μM)K_{R1} = 1(\mu M) Dissociation constant of IPTG and LacI Estimated from [7:2]
n1=2n_1 = 2 Hill coefficient of IPTG and LacI Estimated from [7:3]
r1=4.4E3(s1)r_1 = 4.4E-3(s^{-1}) The degradation rate of mRNAcsgAMBPmRNA·CsgA·MBP Estimated from [8:4]
k1=0.57(s1)k_1 = 0.57(s^{-1}) The translation rate of mRNAcsgAMBPmRNA·CsgA·MBP Estimated from [8:5]
k2=0.44(s1)k_2 = 0.44(s^{-1}) The translation rate of mRNAcsgAMTmRNA·CsgA·MT Estimated from [8:6]
r2=6.3E5(s1)r_2 = 6.3E-5(s^{-1}) The degradation rate of csgAMBPCsgA·MBP Estimated from [9:2]
r3=3.4E3(s1)r_3 = 3.4E-3(s^{-1}) The degradation rate of mRNAcsgAMTmRNA·CsgA·MT Estimated from [8:7]
r4=4.88E5(s1)r_4 = 4.88E-5(s^{-1}) The degradation rate of csgAMTCsgA·MT Estimated from [9:3]
k3=1.8E19(Ls1)k_3 = 1.8E-19(L·s^{-1}) The export rate of CsgA·MBP/MT Estimated from Team:TU_Delft (2015)
KD1=6.655(μM)K_{D1} = 6.655(\mu M) Binding affinity between CsgA·MBP and AgNP(Ag+Ag^+) Estimated from experiment
n2=0.5086n_2 = 0.5086 Hill coefficient between CsgA·MBP and AgNP(silver ions) Estimated from experiment
KD2=0.3559(μM)K_{D2} = 0.3559(\mu M) Binding affinity between CsgA·MT and Cu2+Cu^{2+} Estimated from experiment
KD3=0.1676(μM)K_{D3} = 0.1676(\mu M) Binding affinity between CsgA·MT and Cd2+Cd^{2+} Estimated from experiment
n3=0.4614n_3 = 0.4614 Hill coefficient in the binding of between Cu2+Cu^{2+} and CsgA·MT Estimated from experiment
n4=0.3182n_4 = 0.3182 Hill coefficient in the binding of between Cd2+Cd^{2+} and CsgA·MT Estimated from experiment

  1. 1.  Eberhardt, J., Santos-Martins, D., Tillack, A.F., Forli, S. (2021). AutoDock Vina 1.2.0: New Docking Methods, Expanded Force Field, and Python Bindings. Journal of Chemical Information and Modeling. ↩︎

  2. 2.https://doi.org/10.5281/zenodo.7037337 ↩︎

  3. 3.  H. A. de Boer , L. J. Comstock and M. Vasser, Proc Natl Acad Sci U S A., 1983, DOI:10.1073/pnas.80.1.21 ↩︎

  4. 4.  U. Alon,An Introduction to Systems Biology: Design Principles of Biological Circuits,Chapman and Hall/CRC,Boca Raton,2019 ↩︎ ↩︎ ↩︎

  5. 5.  M. Stamatakis and N. V. Mantzaris, Biophysical journal, 2009, DOI:10.1016/j.bpj.2008.10.028 ↩︎ ↩︎ ↩︎ ↩︎

  6. 6.  T. Kalisky, E. Dekel and U. Alon, Phys Biol., 2007, DOI:10.1088/1478-3975/4/4/001 ↩︎ ↩︎

  7. 7.  S. Semsey, L. Jauffred, Z. Csiszovszki, J. Erdossy, V. Stéger, S. Hansen and S. Krishna, Nucleic acids research, 2013, DOI:10.1093/nar/gkt351 ↩︎ ↩︎ ↩︎ ↩︎

  8. 8.  J. R. Kelly, A. J. Rubin, J. H. Davis, C. M. Ajo-Franklin, J. Cumbers, M. J. Czar, K. de Mora, A. L. Glieberman, D. D. Monie and D. Endy, Journal of biological engineering, 2009, DOI:10.1186/1754-1611-3-4 ↩︎ ↩︎ ↩︎ ↩︎ ↩︎ ↩︎ ↩︎ ↩︎

  9. 9.  M. Halter, A. Tona, K. Bhadriraju, A. L. Plant and J. T. Elliott, Cytometry, 2007, DOI:10.1002/cyto.a.20461 ↩︎ ↩︎ ↩︎ ↩︎

  10. 10.  E. Lee, D. H. Kim, Y. Woo, H. G. Hur and Y. Lim ,Biochem Biophys Res Commun. ,2008, DOI:10.1016/j.bbrc.2008.09.042 ↩︎

  11. 11.  A. L. Lira, R. S. Ferreira, R. Torquato, H. Zhao, M. Oliva, S. A. Hassan, P. Schuck and A. A. Sousa, Nanoscale, 2018, DOI: 10.1039/c7nr06810g ↩︎

  12. 12.  X. Li, Z. Ren, M. J. C. Crabbe, L. Wang and W. Ma, Ecotoxicol Environ Saf. ,2021,DOI:10.1016/j.ecoenv.2021.112512 ↩︎

  13. 13.  A. H. Mohammad, I. Riku, etc (2021) An experiment-based model quantifying antimicrobial activity of silver nanoparticles on Escherichia coli†, RSC Adv., 2017, 7, 56173, DOI: 10.1039/c7ra10495b ↩︎ ↩︎ ↩︎ ↩︎ ↩︎ ↩︎

  14. 14.  Krishnamurthi VR, Niyonshuti II, Chen J, Wang Y (2021) A new analysis method for evaluating bacterial growth with microplate readers. PLOS ONE 16(1): e0245205 ↩︎ ↩︎ ↩︎