Team:UCL/Model

Overview

Our modelling contains two models, flux balance analysis (FBA) and cellular automata (CA), to simulate bacterial plastic degradation and bioelectricity production for desalinating saline water. The FBA models the biomass growth as well as product secretion of our PET-degrading bacteria, Escherichia coli and Pseudomonas putida. Lactate secretion generated by FBA is used in the CA to model the biofilm growth of the exoelectrogen, Shewanella oneidensis (S. oneidensis) and to simulate current generation. Both models will provide insights for synthetic biology strategies to improve the microbial desalination cell (MDC) performance.

Model page Overview diagram

Part 1: Flux Balance Analysis

Highlights of FBA diagram

Flux Balance Analysis (FBA) comprises the first section of our modeling. It is a mathematical modeling technique used to determine the distribution of fluxes across metabolic reactions and pathways in organisms growing at steady state. By calculating the flow of metabolites through the metabolic network of an organism, FBA can help to predict its growth rate, the rate of production of a particular metabolite as well as help guide metabolic engineering strategies (such as single gene deletion) towards maximizing productivity [1].

Our FBA group aimed at simulating and optimizing bacterial growth and product secretion of our plastic degrading species - E. coli BL21 DE3 and P. putida KT2440. We retrieved the respective Genome-Scale Metabolic models (GEMs) from EMBL-EBI, imported them into MATLAB, and performed FBA using the COBRA toolbox. Our work includes three stages:

Fig. 1: 3 stages of FBA diagram
Figure 1. Three stages of our FBA work. Stage one aimed at choosing well curated Genome-Scale Metabolic Models (GEMs) for and P. putida, simulating their growth on preferred and target carbon sources under different oxygen conditions. Stage two concentrated on product secretion and optimization. Stage three centered on producing representative schematics of work from Stage one and two.
Fig. 2: Pathway diagram
Figure 2. Pathway schematic representation of a model of enzymatic plastic degradation coupling microbial electricity generation. The engineered E. coli secretes PET degrading enzymes, PETase and MHETase, to degrade polyethylene terephthalate (PET) into mono-terephthalic acid (MHET) and ethylene glycol (EG) [2][3]. EG enters E. coli via transporters to be further metabolized in the native secondary pathway expressed an operon in E. coli. Metabolites in this pathway are: glycolaldehyde (GA), glycolate (GC), glycoxalate (GLA), tartronate semialdehyde (TS), and glycerate (GL). While GL enters glycolysis and gets converted into pyruvate (Pyr) to support E. coli’s biomass growth, TPA enters P. putida via TPA transporter encoded by genes from Comamonas sp. strain E6 [4]. Our engineered P. putida expresses transgenic genes to degrade TPA into 1,2-dihydroxy-3,5-cyclohexadiene-1,4-dicarboxylate (DCD), next into 3,4-dihydroxybenzoate (PCA) [5]. PCA enters P. putida’ s endogenous β-ketoadipate pathway. Metabolites in this pathway are: β-carboxy-cis,cis-muconate (β-CM), γ-carboxymuconolactone (γ-CML), β-ketoadipate enol-lactone (β-KG EL), β -ketoadipate (β-KG), β -ketoadipyl-CoA (β-KGCoA), succinyl-CoA (SucCoA), acetyl-CoA (AcCoA) [6]. AcCoA enters tricarboxylic acid cycle (TCA) of P. putida, generating lactate for exoelectrogen S. oneidensis. Native pathway in S. oneidensis metabolizes lactate to support biomass growth and electron production [7][8].

Table 1 summarizes the GEM [9][10], single carbon source, and target product of E. coli and P. putida.

Table 1

Stage 1

We focused on investigating and modeling the biomass growth of E. coli and P. putida. Here are the steps we followed:

Research to find the most recent, consensus GEM for E. coli and P. putida and attempted to reproduce biomass growth rates on glucose under aerobic conditions reported in the original publications. This ensured that we had the properly implemented models to begin with.

Analyses:

As presented in Table 2, prediction of E. coli’ s biomass growth is 0.7847/h while Kip et al has reported 0.6/h under the same glucose and oxygen uptake rates [11]. This discrepancy could be explained by the additional thermodynamic constraints applied and the strain difference as E. coli K12 was used in the publication.

he biomass growth simulation of P. putida is consistent with the reported data from Nogales et al, suggesting that the iJN1462 was well curated [10].

Table 2

To further verify, we plotted assumed and predicted growth rate against glucose uptake rate for both species. According to Fig. 3 and 4, values for growth rate assumptions are in proximity to predictions before maxima growth rates were achieved, once again, indicating that our chosen GEMs are well curated.

Fig. 3: <i>E. coli</i> growth rate on glucose
Figure 3. Comparison of assumed and predicted E. coli growth rate against glucose uptake rate. Simulation performed at oxygen uptake rate = -13.1 mmol•gDCW•h, value reported in the publication of E. coli’s GEM.
Fig. 4: <i>P. putida</i> growth rate on glucose
Figure 4. Comparison of assumed and predicted P. putida growth rate against glucose uptake rate. Simulation performed at oxygen uptake rate = -13.5 mmol•gDCW•h, value reported in the publication of P. putida’ s GEM.

Introduce exogenous reactions to E. coli and P. putida for them to carry the PET and TPA degradation mechanism (Fig. 2), respectively. For the convenience of programing, code the PET degradation pathway as intracellular which does not affect our FBA results. Use one repeat for PET. Simulate biomass growth of E. coli on PET and P. putida on TPA with –20 mmol•gDCW•h oxygen uptake rate.

Analyses:

As presented in Table 3, both E. coli and P. putida achieved biomass on our designated carbon sources under aerobic conditions (oxygen uptake rate = -20 mmol•gDCW•h). E. coli grew on PET as the sole carbon source achieving a theoretical growth rate of 0.6543/h while P. putida grown solely on TPA achieved a theoretical growth rate of 0.4468/h. In growth simulations for both organisms, oxygen was found to be the limiting factor. Simulations indicate that both E. coli and P. putida are unable to grow on PET and TPA anaerobically (Fig. 5). As a result, we have changed our coculture design as well as Microbial Desalination Cell (MDC) operation and configuration design accordingly. Please visit Engineering Page for more information.

Table 3
Fig. 5: oxygen conditions
Figure 5. Comparison of anaerobic and aerobic coculture of E. coli and P. putida.

For E. coli, biomass growth is dependent on PET and oxygen uptake rate. As presented in Fig. 6, E. coli’s biomass growth levels off at 0.6543/h. This result is consistent with the 3D plot in Fig. 7, showing that when oxygen uptake rate reaches -20 mmol•gDCW•h, the maximum PET uptake allowed is -20 mmol•gDCW•h, leading to a maximum biomass growth of 0.6543/h. We learned from literature that despite being a facultative anaerobe, E. coli requires sufficient level of oxygen to grow on the degradation product of PET – Ethylene glycol, which enters E. coli’s secondary metabolism and converts into Pyruvate to support growth via the TCA cycle.

Fig. 6: <i>E. coli</i> growth on PET (2D)
Figure 6. E. coli biomass growth simulation on PET. Oxygen uptake rate = -20 mmol•gDCW•h.
Fig. 7: <i>E. coli</i> growth on PET (3D)
Figure 7. E. coli biomass growth simulation on varying PET and oxygen uptake rates.

For P. putida, an obligate aerobe, biomass growth is dependent on TPA and oxygen uptake rate. Fig. 8 compares P. putida’ s biomass growth on TPA and glucose, the default carbon source in GEM. For this simulation, the oxygen uptake rate was set by default to -13.5 mmol•gDCW•h. The maximum growth on TPA, 0.3/h, is significantly lower than that on glucose, 0.71/h since TPA is not P. putida’ s preferred carbon source. Upon discussion, we decided to pursue a theoretical investigation of P. putida growth on TPA at higher oxygen uptake rates. Therefore, we increased oxygen uptake rate to -20 mmol•gDCW•h and performed growth simulation again. The maximum TPA uptake rate and growth rate are -6 mmol•gDCW•h and 0.4468/h, respectively (Fig. 9).

Fig. 8: <i>P. putida</i> growth on TPA (2D)
Figure 8. Comparison of P. putida growth rate on glucose and TPA. Oxygen uptake rate = -20 mmol•gDCW•h.
Fig. 8: <i>P. putida</i> growth on TPA (3D)
Figure 9. P. putida biomass growth simulation on varying TPA and oxygen uptake rates.

Stage 2

We focused on simulating and optimizing the TPA secretion rate of E. coli and the lactate secretion rate of P. putida. Here are the steps we followed:

Simulate TPA and lactate secretion rates at -20 mmol•gDCW•h oxygen uptake rate. Explore the correlation amongst production rate, biomass growth rate, and carbon source uptake rate.

Analyses:

Overall, E. coli’s production of TPA is sufficient for P. putida’ s growth and lactate secretion. The second section of our team’s modeling – Cellular Automata (CA) has proved that simulated maximum lactate secretion rate was sufficient for the biofilm formation of S. oneidensis. Therefore, we have validated the feasibility of coupling enzymatic plastic degradation to microbial electricity generation.

It has been discovered that TPA secretion of E. coli is proportional to PET uptake rate, however, disproportional to biomass growth although a minimum growth rate of 0.2/h is suggested to base our simulations on in order to correlate to reality. Our 3D plot in Fig. 10 clearly demonstrates that one repeating unit of PET polymer (C10H8O4) is converted into one TPA molecule (C8H6O4).

Fig. 10: <i>E. coli</i> TPA production (3D)
Figure 10. E. coli TPA production simulation on varying PET uptake rate and biomass growth rate.

Simulations on P. putida reveal that lactate secretion rate is inversely related to biomass growth while directly related to TPA uptake rate (Fig. 11, 12). Lactate secretion rate reaches its theoretical maximum 8.2113 mmol•gDCW•h when growth rate is 0.2/h and TPA uptake rate is -7.5 mmol•gDCW•h (Fig. 12).

Fig. 11: <i>P. putida</i> lactate secretion rate vs. biomass growth rate
Figure 11. P. putida lactate secretion simulation 3D plot from the projection of lactate secretion rate against biomass growth only.
Fig. 12: <i>P. putida lactate simulation (3D)</i>
Figure 12. P. putida lactate secretion simulation on varying TPA uptake rate and biomass growth rate.

From analyzing Fig. 12 and receiving advice from our FBA supervisors, we chose the simulate product secretion rate of both E. coli and P. putida at 0.2/h growth rate which is suggested to be the minimum growth for our model to still correlate reality whilst maximizing production secretion.

Perform single gene deletion on the P. putida model to identify, for each organism, a single gene that could maximize product secretion if knocked out. Not necessary to do this for E. coli since TPA does not come from E. coli’s native metabolism.

Analyses:

It has been identified that when gene 836 in P. putida’ s GEM was knocked out, the maximum lactate secretion rate reached 8.2132 mmol•gDCW•h. We learned that this gene is involved in four reactions and encodes iron-sulfur cluster assembly protein. Future work could gear towards understanding the impact of this gene knockout on reaction fluxes in the P. putida model.

Table 4 summarizes the key product secretion results from our modeling at Stage Two.

Table 4

Stage 3

We were mainly working on generating and analyzing 2D and 3D plots as well as drawing pathway schematic. Please visit section “FBA Advice for Future Teams” on our Contribution Page to check out resources we have used for graphics.

Part 2: Cellular Automata

A Cellular Automata (CA) is a grid of a specific size where each node has a defined rule set based on cell behaviours. From a relatively simple rule set we can model cell interactions and predict biofilm behaviours. The cells on the grid will have one of a finite number of states that will be updated over discrete time steps according to rules based on the state of the adjacent neighbouring cells.

Model Aims:

  • To investigate growth and behaviour of S. oneidensis biofilm on the anode for electricity generation.
  • To predict desalination efficiency achieved over two different current output models.
  • To inform synthetic biology strategies for engineering S. oneidensis to improve electron transfer and maximise desalination efficiency.

Overview

We developed a CA model of a 3D grid on the anode surface with each node representing a cluster of S. oneidensis cells in the biofilm which can have one of 4 states: empty, alive, quiescent or dead. The state of each cell in the grid is updated by rules based on the substrate concentration and state of the neighbouring cells. The diffusion of the substrate and the electron shuttling molecule (A flavin, referred to as the mediator) in the 3D grid were modelled and used to compute the equations for updating the CA. The resulting grid of cells was used to model the current generated by extracellular electron transfer and linked to a desalination efficiency model.

Cellular automata model overview diagram
Figure 13. Overview of the 4 stages in cellular automata model. Stage 1: Models the diffusion of the substrate, lactate (blue) and the oxidized (orange) and reduced forms (green) of the mediator in the 3D grid. The lactate concentrations at the boundary is determined by the lactate production of the co-culture modelled by the FBA. Stage 2: The cellular automata which is updated based on lactate availability and the number of empty neighbours of each cell. Stage 3: The cellular automata is linked to an extracellular electron transfer model of the biofilm to generate a current output. Stage 4: A desalination model that converts the current generated in stage 3 into a desalination rate.

Stage 1: Substrate Diffusion Model

The first part of the Cellular automata is the diffusion of the substrates which is used to determine the concentrations of lactate and the oxidized and reduced forms of the mediator in the grid over time.

Diffusion was modelled with Fick’s second law in 3 dimensions using a central finite difference approximation which calculates the substrate concentration at each point in the grid using the concentrations of the six neighbouring points. The grid is also updated in time using the Forward Euler method.

\[\frac{{\partial C}}{{\partial x}} = D\left( {\frac{{{\partial ^2}C}}{{\partial {x^2}}} + \frac{{{\partial ^2}C}}{{\partial {y^2}}} + \frac{{{\partial ^2}C}}{{\partial {z^2}}}} \right)\]
  • C = concentration of substrate (mM)
  • D = Diffusion coefficient (m2s-1)
  • x, y, z = spatial dimensions of 3D grid

The finite difference approximation of the equation using Forward Euler is

\[\frac{{C_{x,y,z}^{t + 1} - C_{x,y,z}^t}}{{\Delta t}} = D\left( {\frac{{C_{x + 1,y,z}^t - C_{x - 1,y,z}^t + C_{x,y + 1,z}^t - C_{x,y - 1,z}^t + C_{x,y,z + 1}^t - C_{x,y,z - 1}^t - 6C_{x,y,z}^t}}{{{h^2}}}} \right)\]

Assumptions:

  • Diffusion is 50% slower in the biofilm than the bulk liquid as biofilm is denser
  • The chamber is well mixed so there are no concentration gradients within the bulk liquid

As mentioned above, lactate diffusion updates are computed through solving the heat equation using the Forward Euler method. The Forward Euler method uses the gradient, which signals how rapidly a variable is changing over time, to estimate the value of the variable based on a given starting point. The key difference is that a differential equation depicts a variable’s rate of change, but not the value of the variable itself. Therefore, the Forward Euler is an approximation of the actual value based on its gradient over time. This approximation becomes more accurate as the time-steps between one estimate to the next get shorter. The heat equation is a partial differential equation that can be used for simulating the way particles diffuse in a fluid, and in our context how lactate diffuses in the anode chamber.

However, the shorter these time steps are, the more computation is needed to run the simulation. In order to find a reasonable compromise between accuracy and computational time, we can use Von Neumann Stability Analysis. This method can find what is the time-step size required so that errors from one step to the next does not grow non-linearly, as in, add up together exponentially, leading to inaccurate results. Below is the resulting equation depicting tilme-step size for the heat equation in 3 dimensions, where:

\[\Delta t \le \frac{{6D}}{{{h^2}}}\]
  • h represents the length of each side cube in the simulated grid, in metres. This is estimated by how many cells we are simulating per cube, and will be discussed below.
  • D represents the diffusion coefficient of lactate in water, with units of m2s-1.

Each S. oneidensis cell is about 3µm in length (retrieved from Bionumbers). In our simulations we estimated each cell to occupy a cube with a side length of 5µm, as we are assuming some extra space is taken up by extracellular matrix. If we aimed to simulate each cell in the model as an agent, h would equal 5µm. Based on the stability analysis shown above, this resulted in a timestep of approximately 0.002 seconds. If we intend to simulate a week of cell growth, more than 300 million lactate diffusion updates must be computed, which leads to an unrealistic time required for computation.

In order to circumvent this, we can have each agent in the grid represent a cluster of cells. Through this, we can increase the value of h, which drastically lengthens our time-steps. However, using each agent to represent a cluster of cells is not optimal, since a true agent based model would aim to simulate each cell individually. The compromise that we reached between these two factors was having each cube in the grid have a side length of 25µm, 5 times the original value of h, and therefore representing 125 cells per agent. This new value of h renders the new time-step to be approximately 0.06 seconds, leading to 12 million lactate diffusion updates per week, which becomes a manageable task.

We have modelled a representative area of the anode of finite size L x L x L where L is the number of nodes in each plane of the grid. This is a modelled part of a larger system and at the boundaries of the grid, we need to account for this, so they have been defined to be periodic for all species modelled which means that the flux at the boundaries are equal.

\[\frac{{\partial C_{x = 0,y,z}^t}}{{\partial x}} = \frac{{\partial C_{x = L,y,z}^t}}{{\partial x}}\] \[\frac{{\partial C_{x,y = 0,z}^t}}{{\partial y}} = \frac{{\partial C_{x,y = L,z}^t}}{{\partial y}}\]

The boundary at the anode surface is defined as an impermeable surface (insulated boundary) as there will be no flow of substrate through it.

\[\frac{{\partial C_{x,y,z = 0}^t}}{{\partial z}} = 0\]

Lactate Diffusion

The initial lactate concentration at t = 0 and boundary condition of the bulk liquid at the top interface of the grid (z = L) is defined by the output of the FBA.

Given by the FBA, the plastic degrading coculture could generate 8.2132 mmol lactate/hour/grams of Ash Dry Cell Weight (ADCW). As explained in the FBA section, some of the carbon obtained from PET has to go towards biomass growth, and some towards lactate production. The current steady state simulated by the FBA for each bacterium in the coculture is of 0.2 grams. This sums up to 0.4 grams of coculture.

Based on this, we can calculate how many millimoles of lactate flow into the anode chamber per unit time. The units of time are based on how long (seconds) each timestep is regarding the lactate diffusion update. The length of this timestep is calculated in each simulation. Its based on the Von Neumann stability analysis (See above), which is influenced by the size of each cell in the grid and the speed at which lactate diffuses (Diffusion coefficient of lactate). Knowing the lactate influx, we can then calculate how much lactate reaches the grid at the surface of the anode that we are simulating.

As explained above, the toroidal boundary conditions aim to allow us to compute biofilm growth and current density through extrapolating a grid of size L over the whole anode surface. Given this size, we can calculate how much volume the simulated anode surface takes up compared to the volume of the chamber. We will call "bulk liquid" the volume of the chamber that is not part of the CA simulation. By assuming that lactate is evenly distributed throughout the bulk liquid, we can use the concentration in the bulk liquid as a boundary condition for the lactate diffusion grid at the start of each update. This will generate a concentration gradient towards the anode surface.

How much lactate flows onto the grid from the bulk liquid is then built back into the overall lactate concentration in the rest of the chamber. This process is repeated during every lactate update, computing the balance between lactate flowing into the chamber and being consumed by S. oneidensis .

The steady-state initial and boundary conditions can be stated generally as equations, f(t) and g(t) that are function of time.

Initial condition at t = 0:

\[C_{x,y,z}^{t = 0} = f(t)\]

Boundary condition at z = L

\[C_{x,y,z = L}^t = g(t)\]

Mediator Diffusion

Initial condition at t = 0:

$$C_{x,y,z}^{t = 0}({M_{ox}}) = 0$$ $$C_{x,y,z}^{t = 0}({M_{red}}) = 1mM$$

Boundary condition at z = L:

$$C_{x,y,z = L}^t({M_{ox}}) = 1mM$$ $$C_{x,y,z = L}^t({M_{ox}}) = 1mM$$
  • Mox = oxidized mediator
  • Mred = reduced mediator
Summary of diffusion boundary conditions
Figure 14. Summary of the boundary conditions of lactate and the oxidized and reduced mediator in the 3D grid. The anode (dark grey) at z = 0 has no flux at the boundary as defined by the equation given. The grey and clear surfaces represent the periodic boundaries where the flux at one side is equal to the flux on the other side.

Stage 2: Cellular automata biofilm growth model

The substrate consumption was modelled using double-limited Monod kinetics which assumes that the lactate uptake rate is controlled by the availability of lactate and the electron acceptor [12,13]. Since the model consists of two modes of electron transfer: direct transfer, and mediated transfer. The electron acceptor will either be the anode or the mediator. The mediator is transfer limiting as it is present at low concentrations and in both redox states in the biofilm. Meanwhile the redox potential between the terminal donator and the anode controls the rate of electron generation at the anode [12].

Therefore, two rates of electron generation are present in the biofilm, those in direct contact with the anode and, those with a mediated contact. Two lactate substrate consumption equations were used depending on the position of the cells in the biofilm. The consumption rate was computed from the specific growth rate of each cell.

For the cells in the first layer of the grid in direct contact with the anode, the substrate consumption rate will depend on the cell outer membrane potential:

\[q = \frac{{{\mu _{\max }}X}}{Y}\left( {\frac{{{C_{lactate}}}}{{{C_{lactate}} + {K_{lactate}}}}} \right)\left( {\frac{1}{{1 + \exp \left[ {\frac{{ - F}}{{RT}}\left( {E - {E_{KA}}} \right)} \right]}}} \right)\]

The substrate consumption rate of the cells not in direct contact is defined by the availability of the oxidized form of the mediator which will be accepting the electrons generated.

\[q = \frac{{{\mu _{\max }}X}}{Y}\left( {\frac{{{C_{lactate}}}}{{{C_{lactate}} + {K_{lactate}}}}} \right)\left( {\frac{{{C_{{M_{ox}}}}}}{{{C_{{M_{ox}}}} + {K_{{M_{ox}}}}}}} \right)\]

where

  • q = substrate utilization rate (mmol -1)
  • µmax = maximum growth rate (s-1)
  • X = cell density (gm-3)
  • Y = yield of biomass on lactate (g DW mol-1)
  • Klactate = lactate half saturation constant (mM)
  • KMox = oxidized mediator half saturation constant (mM)
  • EOM = cytochrome (outer membrane) potential (V)
  • EKA = potential at half the maximum current density (V)
  • F = Faraday’s constant (C mol-1)
  • R = ideal gas constant (J mol-1 K-1)
  • T = temperature (K)

We adapted Conway’s Game of Life [14], a CA that evolves according to the number of live neighbours of each cell, to develop rules for updating the grid based on mathematical equations. Each cell can have one of 4 states which is updated in each time step using the results calculated from the previous stages.

State of cells in 3D grid: live, quiescent, dead or empty
Figure 15. The 4 states of cell in the CA grid, with the relevant properties of each state.

Assumption: Growth is limited by the availability of the substrate and the electron acceptor (anode or mediator depending on electron transfer mechanism)

Growth update rules:

  1. Live cells can only divide into an empty neighbouring cell which is selected randomly
  2. A cell becomes quiescent if there is below a minimum amount of lactate present
  3. Quiescent cells can become alive if they receive above a certain amount of lactate
  4. A cell can die after staying quiescent for 24 hours

Growth Update steps:

  1. Calculate the growth rate of alive and quiescent cells in grid using the substrate uptake rate
  2. Calculate the probability of division using the growth rate
  3. Compute the probability of alive cells becoming quiescent
  4. Compute probability of quiescent cells dying
  5. Update cells in grid using growth update rules
State of cells in 3D grid: live, quiescent, dead or empty
Figure 16. Decision flow diagram of the rules of the CA and how the grid is updated in each time step.

Stage 3: Current Output Model

The objective of the current output model is to simulate electricity generation by the exoelectrogen, S. oneidensis by building on the cellular automata to model electron transfer from the biofilm to the anode. Direct and mediated electron transfer have been modelled as these mechanisms have observed in S. oneidensis.

Mechanisms of electron transfer by S. oneidensis
Figure 17. The different mechanisms of extracellular electron transfer (EET) mechanisms. a. direct electron transfer that occurs as a result of close contact between the cell's outer membrane proteins (referred to as cytochromes) and the anode; b. mediator electron transfer – oxidised flavin mononucleotide mediator is reduced by the cell and transports electrons across in a shuttle manner. Mr = mediator in reduced state; Mo = mediator in oxidised state

Model assumptions

  • All electrons are transferred to the anode
  • The rate limiting step determines current output as reactions occur in series
  • There is a constant potential loss in the biofilm over time

The layer of the biofilm in direct contact with the electrode surface transfer electrons from the cell cytochromes to the anode. As the electron don’t have to travel over a long distance in this mechanism, direct transfer has the lowest extracellular potential loss [15]. This mechanism contributes less than 25% of the current generated by S. oneidensis biofilms, therefore in our model the effect of the potential losses on the total current output were assumed to be negligible [16].

Assumptions:

  • Only the cells in the first layer of the grid in contact with the anode undergo direct electron transfer
  • Only live cells contribute to the electron generation
  • Negligible potential loss between the electrode interface and anode so current generated is close to the maximum current density

Maximum Current Density

The current generated by the direct contact mechanism is calculated using the equation for the maximum current density generated by the monolayer. The maximum current density describes the highest possible current that can be generated by the biofilm without any potential losses [15].

\[{j_{\max }} = {\gamma _s}{q_{\max }}{X_f}{L_f}\]

Where:

  • γs = conversion factor from mass of substrate to coulombs (C)
  • qmax = maximum specific rate of substrate utilization (mol g−1 dwt s−1)
  • Xf = concentration of active biomass in biofilm (g dwt m−3)
  • Lf = biofilm thickness active in electron transfer (m)

Cells further away from the anode surface transfer their electrons over a long distance through two main mechanisms: conductive pili/nanowires and soluble electron shuttles. This model focuses on the latter as it makes up 75% of the electricity generated by S. oneidensis [16].

The equation of the redox reaction driving the mediator electron transfer mechanism is:

\[{M_{ox}} + 2e \mathbin{\lower.3ex\hbox{$\buildrel\textstyle\rightarrow\over {\smash{\leftarrow}\vphantom{_{\vbox to.5ex{\vss}}}}$}} {M_{red}}\]

The processes of this mechanism included in the model are (i) secretion of the mediator by the cells undergoing MET, (ii) Electron transfer to the mediator by the biofilm, (iii) Electron generation at the anode through reversible reactions and (iv) the current density from the most limiting process of electron generation at anode surface and long-range diffusion.

Mechanism of mediator transfer
Figure 18. Mechanism of mediator electron transfer by exoelectrogen. The cells in the biofilm secrete electrons which are transferred to the oxidized mediator generating the reduced form (1). The mediator then diffuses through the biofilm to reach the anode surface (2) where the electrons are delivered to the anode in a reversible reaction that produces the oxidized mediator (3).

Assumptions:

  • Oxidized and reduced mediator have the same diffusion coefficient and standard heterogeneous rate constant
  • Mediator movement is only affected by diffusion biofilm migration or electrochemical interactions).
  • The rate of mediator secretion depends on the lactate concentration
  • Mediator is only secreted by live cells undergoing MET
  • Electron transfer from the mediator at anode surface is reversible
  • Electron transfer from the biofilm to the mediator is irreversible

Mediator Secretion

\[{R_s} = {S_{\max }}\left( {\frac{{{C_{lactate}}}}{{{C_{lactate}} + {K_{lactate}}}}} \right)\]
  • Rs= specific mediator secretion rate (mmol cell-1 s-1)
  • Smax = maximum mediator secretion rate (mmol gDW-1 h-1)
  • Clactate = Lactate concentration (mM)
  • Klactate = Lactate half saturation constant (mM)

Electron Transfer by Biofilm

The model of the electron transfer in the biofilm was adapted from Renslow at al. [12]. The concentration of the mediator in both the reduced and oxidized state was evaluated for each cell active in mediator electron transfer (cells above the first layer). The transfer rate is described in terms of the fraction of electrons transferred by each mediator and the lactate consumption rate of each cell.

\[{R_M} = \frac{y}{n}fq\]
  • RM = reduction rate of oxidized mediator in biofilm (mol m−3 s−1)
  • y = electron equivalence of the substrate (unitless)
  • f = fraction of electrons recoverable for current
  • n = number of electrons transferred by a redox mediator
  • q = specific substrate utilization rate via MET (mol m−3 s−1)

Electron Generation at the Anode

As the reaction at the anode surface is reversible, the net rate of electron generation is determined by the difference between the rate of the forward and reverse reaction. With the anode being the electron acceptor, the rate is determined by the anode potential as well as the reaction on the anode surface.

\[{r_M} = {k_0}\left( {C{{_{{M_{red}}}^t}_{x,y,z = 0}}\exp \left[ {(1 - \alpha )\frac{{nF}}{{RT}}\left( {E - {E^0}} \right)} \right] - C{{_{{M_{ox}}}^t}_{x,y,z = 0}}\exp \left[ { - \alpha \frac{{nF}}{{RT}}\left( {E - {E^0}} \right)} \right]} \right)\]
  • rM = rate of electron generation at the anode (mol m-2 s-1)
  • k0standard heterogeneous rate constant for the mediator redox reaction (ms-1)
  • α = the transfer coefficient for the mediator redox reaction
  • E0 = the standard redox potential for the mediator redox reaction (V)
  • E = anode potential (V)

Computation of mediator concentration

  • The mediator is reduced by the biofilm and oxidized at the electrode at each iteration of the CA
  • The rate of change in the concentrations of the oxidized and reduced mediator is determined by the mass balance of the sum of the rates calculates.

Current Density Output

As the reactions happen in series, the current density will be limited by the slowest process. Over time, as the biofilm gets thicker, it is expected that the diffusion of the mediator through the cells in the biofilm will become limiting however, the electron generation at the anode surface is initially the limiting process. Therefore, the current output is defined as

\[j = nF{r_{\lim }}\] \[{r_{\lim }} = \frac{{{D_M}\Delta {C_M}}}{{\Delta z}},{r_{\lim }} = {r_M}\]
  • rlim = limiting reaction
  • DM = mediator diffusion coefficient (m2s-1)
  • ΔCM= concentration gradient of mediator (mM)
  • Δz = diffusion distance (m)

Stage 4: From current output to desalination rate

The Nominal desalination rate (L m2 h-1) represents the volume of freshwater per square meter of cell membrane per desalination cycle and includes normalised current output [17]:

\[NDR = \frac{{\frac{1}{{SEP}}{E^0}\int {j(t)dt} }}{t}\]

Assumptions: SEP estimated from max current density ever produced (8.4 A m2) and max SEP (0.4 kWh m-3) [18], and comparing to own current density produced.

Results and conclusions

The aim of the simulations was to find the right initial conditions towards reaching a stable steady state for S. oneidensis growth, as well as how much desalination this steady state could carry out. We considered a stable steady state when the number of active cells did not drastically fluctuate and this behaviour was maintained for at least one week. The size of this steady state is determined by the rate of lactate influx generated by the plastic-degrading co-culture, as well as the rate of lactate consumption by S. oneidensis . An essential factor for reaching this steady state was the initial concentration of lactate in the chamber as S. oneidensis was introduced. This depended on how long the co-culture had been degrading plastic for and how much lactate flowed into the anode chamber. We will refer to this time length as accumulation time. If the accumulation time was too low, the cells would not have enough lactate to enter exponential growth and will die before reaching a steady state. On the other hand, if the accumulation time was too high, the cells would grow unsustainably, consume all of the lactate in the chamber, and enter an aggressive death phase that could not be rescued by the lactate influx.

As mentioned above, the lactate influx depended on the size of the plastic-degrading co-culture and on its ability to generate lactate from PET. The exploration depicted below was based on small co-culture, of 0.4g in ADCW that was able to generate 3.2853 mmol of lactate per hour. One week of accumulation time for this co-culture reliably lead to a stable steady state of S. oneidensis growth. Below is a graph showing Active (Blue), Quiescent (Orange) and Dead cells (green) over time during an initial week of cell growth. Important to note that each agent represents 125 cells.

Results of active, quiescent and dead cells
Figure 19. Agent type over a week, with time in a hours,showing active cells in blue, quiescent cells in orange, and dead cells in green. Each agent behaves as 125 cells.

As observed, the cells start to grow exponentially over the first half of the week. However, the resulting biofilm consumes lactate rapidly, which depletes the chamber. This leads to a rapid increase in the number of quiescent cells, which in turn kickstarts a sigmoidal curve of dead cells over time. The continuous spiking pattern is a result of the tight balance between lactate consumption and lactate influx, causing cells to switch over and back from an active state to a quiescent state under the fluctuations in lactate concentration. This spiking pattern stabilises over time, as the number of active cells approaches a stable steady state and the rates of quiescence decrease.

Below are a series of graphs showing the following week of the simulation, depicting cell behaviour, lactate concentration and different mechanisms of current density. The model has not reached a complete steady state yet, as although the number of active cells is not drastically fluctuating (see figure 2a), the lactate concentration is still decreasing over time (see figure 2f). The expected behaviour from this point is that the ratio of quiescent cells will increase over time, leading to a small bout of cell death that will result in the stable steady state.

A
B
C
D
E
F
G
H
Figure 20. Cell behaviour, lactate concentration and current density results of the second week of simulation. A) Agent type over time B) Total current density over time C) Mediator concentration over time D) Maximum direct transfer current density over time E) Biofilm thickness over time F) Lactate concentration over time G) Oxidised vs Reduced mediator over time H) Mediator transfer current density over time

As explained above, we can convert current density into desalination rate, with units of litres per metre squared per hour. In order to do so, we must integrate our current density graph, and therefore we need to carry out a polynomial fit. A second degree polynomial fit to our current density graph is shown below.

Current density over time
Figure 21. Current density over time, shown above, compared to polynomial approximation, shown below.

After integrating the second degree polynomial shown above, we arrive at an average rate of 0.623 litres per metre squared per hour of desalinated water. This value is reached in a 250 ml anode chamber, with a 0.4g co-culture. This poses interesting questions regarding scale-up. It is also of interest to note that no direct transfer-based current was generated, due to all cells dying at the anode surface. A higher lactate influx rate is likely to prevent this from happening, as the biofilm will not have to withstand an aggressive death phase. This model will be used to explored the scenarios described above, and guide next year’s UCL iGEM team in developing a real-life Microbial desalination cell in the laboratory.

Table of parameters used in model

Databases, Tools, Software

  1. EMBL-EBI BioModels
  2. MATLAB
  3. COBRA Toolbox
  4. PyCharm CE
  5. BioNumbers

Graphic Design

  1. Biorender
  2. Pixton.com
  3. Draw.io
  4. Lucidcharts
  5. Powerpoint

Literature

  1. Orth JD, Thiele I, Palsson BØ. What is flux balance analysis? Nat Biotechnol [Internet]. 2010 Mar;28(3):245–8. Available from: https://pubmed.ncbi.nlm.nih.gov/20212490
  2. Yoshida S, Hiraga K, Takehana T, Taniguchi I, Yamaji H, Maeda Y, et al. A bacterium that degrades and assimilates poly(ethylene terephthalate). Science (80- ) [Internet]. 2016 Mar 11;351(6278):1196 LP – 1199. Available from: http://science.sciencemag.org/content/351/6278/1196.abstract
  3. Salvador M, Abdulmutalib U, Gonzalez J, Kim J, Smith AA, Faulon J-L, et al. Microbial Genes for a Circular and Sustainable Bio-PET Economy. Genes (Basel) [Internet]. 2019 May 16;10(5):373. Available from: https://pubmed.ncbi.nlm.nih.gov/31100963
  4. Hosaka M, Kamimura N, Toribami S, Mori K, Kasai D, Fukuda M, et al. Novel tripartite aromatic acid transporter essential for terephthalate uptake in Comamonas sp. strain E6. Appl Environ Microbiol [Internet]. 2013/08/02. 2013 Oct;79(19):6148–55. Available from: https://pubmed.ncbi.nlm.nih.gov/23913423
  5. Johnson CW, Beckham GT. Aromatic catabolic pathway selection for optimal production of pyruvate and lactate from lignin. Metab Eng [Internet]. 2015;28:240–7. Available from: http://www.sciencedirect.com/science/article/pii/S1096717615000075
  6. Wang J-Y, Zhou L, Chen B, Sun S, Zhang W, Li M, et al. A functional 4-hydroxybenzoate degradation pathway in the phytopathogen Xanthomonas campestris is required for full pathogenicity. Sci Rep [Internet]. 2015 Dec;5(1):18456. Available from: http://www.nature.com/articles/srep18456
  7. Kane AL, Brutinel ED, Joo H, Maysonet R, VanDrisse CM, Kotloski NJ, et al. Formate Metabolism in Shewanella oneidensis Generates Proton Motive Force and Prevents Growth without an Electron Acceptor. Silhavy TJ, editor. J Bacteriol [Internet]. 2016 Apr;198(8):1337–46. Available from: https://jb.asm.org/content/198/8/1337
  8. Luo S, Guo W, H. Nealson K, Feng X, He Z. 13C Pathway Analysis for the Role of Formate in Electricity Generation by Shewanella Oneidensis MR-1 Using Lactate in Microbial Fuel Cells. Sci Rep [Internet]. 2016 Aug;6(1):20941. Available from: http://www.nature.com/articles/srep20941
  9. Kim H, Kim S, Yoon SH. Metabolic network reconstruction and phenome analysis of the industrial microbe, Escherichia coli BL21(DE3). Virolle M-J, editor. PLoS One [Internet]. 2018 Sep;13(9):e0204375. Available from: https://dx.plos.org/10.1371/journal.pone.0204375
  10. Nogales J, Mueller J, Gudmundsson S, Canalejo FJ, Duque E, Monk J, et al. High‐quality genome‐scale metabolic modelling of Pseudomonas putida highlights its broad metabolic capabilities. Environ Microbiol [Internet]. 2020 Jan;22(1):255–69. Available from: https://onlinelibrary.wiley.com/doi/abs/10.1111/1462-2920.14843
  11. Kiparissides A, Hatzimanikatis V. Thermodynamics-based Metabolite Sensitivity Analysis in metabolic networks. Metab Eng [Internet]. 2017 Jan;39:117–27. Available from: https://linkinghub.elsevier.com/retrieve/pii/S109671761630218X
  12. Martin E. John Conway’s Game of Life [Internet]. [cited 2020 Oct 24]. Available from: https://bitstorm.org/gameoflife/
  13. Renslow R, Babauta J, Kuprat A, Schenk J, Ivory C, Fredrickson J, et al. Modeling biofilms with dual extracellular electron transfer mechanisms. Phys Chem Chem Phys [Internet]. 2013 Nov 28 [cited 2020 Oct 24];15(44):19262–83. Available from: /pmc/articles/PMC3868370/?report=abstract
  14. Ledezma P, Greenman J, Ieropoulos I. Maximising electricity production by controlling the biofilm specific growth rate in microbial fuel cells. Bioresour Technol. 2012 Aug 1;118:615–8.
  15. Search BioNumbers - The Database of Useful Biological Numbers [Internet]. [cited 2020 Oct 24]. Available from: https://bionumbers.hms.harvard.edu/search.aspx
  16. Picioreanu C, Head IM, Katuri KP, van Loosdrecht MCM, Scott K. A computational model for biofilm-based microbial fuel cells. Water Res. 2007 Jul 1;41(13):2921–40.
  17. Ramírez-Moreno M, Rodenas P, Aliaguilla M, Bosch-Jimenez P, Borràs E, Zamora P, et al. Comparative Performance of Microbial Desalination Cells Using Air Diffusion and Liquid Cathode Reactions: Study of the Salt Removal and Desalination Efficiency. Front Energy Res [Internet]. 2019 Dec 5 [cited 2020 Oct 24];7:135. Available from: https://www.frontiersin.org/article/10.3389/fenrg.2019.00135/full
  18. Kim Y, Logan BE. Microbial desalination cells for energy production and desalination. Desalination. 2013 Jan 2;308:122–30.


Please visit our Github page linked in the footer below for code used in our Modeling!

Our Sponsors
UCL School of Engineering UCL Biochemical Engineering Genscript
© 2020 UCL iGEM