Team:Toronto/Model
iGEM Toronto 2018 Mathematical Modelling
B. Kell, F. Vermehr, Q. Vilchez, P. Malik
October 2018
Contents
1
Introduction
2
2
A Differential Bioreactor Model
2
2.1
Goals
2
2.2
Mathematical Formulation
2
2.2.1
Concepts
2
2.2.2
The Model
3
2.2.3
Bioreactor Assumptions
4
2.3
Equations
4
2.4
Numerical Solution
5
2.5
Sensitivity Analysis
5
2.5.1
The Effect of Dilution Rate on Bioreactor Performance
6
2.5.2
The Effect of Initial Metal Concentration on Bioreactor Performance
7
2.5.3
The Effect of Attachment Rate on Bioreactor Performance
7
2.5.4
The Effect of Detachment Rate on Bioreactor Performance
8
2.5.5
The Effect of Initial Cellular Population on Bioreactor Performance
8
2.5.6
The Effect of Maximum Specific Growth Rate on Bioreactor Performance . .
9
2.6
Hypothetical Applications
9
3
A Mechanistic Buoyancy Model: Maximum Carrying Capacity Determination
10
3.1
Goal
10
3.2
Mathematical Formulation
10
3.2.1
Concepts
10
3.2.2
Assumptions
11
3.2.3
The Model
11
3.3
Analytic Solution
11
3.4
Determining Buoyant Force: Temporal Tracking Algorithm
12
3.5
Discussion On Maximum Carrying Capacity
14
4
Growth Dynamics
15
4.1
Objective
15
4.2
Approach
15
4.3
Results
15
5
Discussion and Future Directions
17
6
Appendix
19
7
References
19
1
1
Introduction
Our goal in the dry lab this year was to create four different models that allows our wet lab team to
characterize their results, and allow future researchers to benchmark their results, creating standard
measures in the field of cellular flotation. First we created a generic differential bioreactor model that
allowed our team to predict the effectiveness of our E. coli cells to clean waste waters if coupled with
any surface binding method. We performed a complete sensitivity analysis on this model to allow
future researchers to reuse this model with completely different parameters, strain of cells and object
of waste. Then we created an algorithm that can track cellular flotation from frame to frame, and
characterize exactly how the cells float; previously, we could only tell whether they floated or not.
This coupled with our ODE buoyancy model allows us to define a maximum carrying capacity for
each strain. Both of these models allowed our team to benchmark their results, and will allow future
researchers to quantify the performance of flotation as well.
Note: In the very end of this paper, we included a nomenclature defining all the variables used.
2
A Differential Bioreactor Model
2.1
Goals
• Explore a possible application for our genetically engineered E. coli biomass that utilizes flota-
tion.
• Develop a generic bioreactor that can be reused in many different conditions and for a variety
of purposes.
2.2
Mathematical Formulation
2.2.1
Concepts
Monod Equation [7]
In order to run a bioreactor, we need to understand at what rate the E. coli grow, and how they
respond to environmental conditions. Bacterial growth can be separated into four phases:
1. The Lag Phase: Little growth is observed.
2. The Exponential Phase: Exponential growth is observed after the cells get used to the environ-
ment, and until the limiting substrate is still in surplus.
3. The Stationary Phase: No growth is observed, death rate and growth rate are approximately
equal.
4. Death Phase: There is no more substrate, thus the population dies.
If N is the concentration of biomass population, then the exponential growth phase is given by:
N
= µN
where µ is the specific growth rate, the ”rate of increase of cell concentrations per unit cell concen-
tration” (h−1) [7].
µ depends on the substrate concentration; the relation is given by the Monod
equation:
S
µ(S) = µm
KS + S
2
where S is the substrate concentration (in mol L−1) , and KS is the Monod Constant (in mol L−1).
As substrate is consumed, more biomass is being created at the rate of:
dN
= −γdS
dS
dt
where t is time. γ should be interpreted as “the ratio of the mass of cells formed to the mass of
substrate consumed” [7].
2.2.2
The Model
Figure 1: The Bioreactor Design
Call the ’Main Reactor’ (M.R) and the ’flotation Tank’ (F.T). Let Xu be the concentration of E. coli
(in mol L−1), Xw the concentration of metal bound to the surface of the E. coli (in mol L−1). Below
we will outline the steps of the bioreactor seen in Figure 1.
1. All the substrate, waste and cellular population is homogeneously distributed between the F.T
and the M.R. Initially there is no metal bound to any of the cells, i.e Xw = 0.
2. When activated, the bioreactor begins to pump solution from the M.R to the F.T at rate D,
this is represented by D(N + Xu + Xw + S).
3. In the F.T, the cells are floating on the surface, and are scraped off and removed. Thus only
Xu and S remain within the solution which is pumped back into M.R at rate D. In practice, a
cellular filter alike the one used for skimming could be placed at the exit of the F.T to capture
any dead cells.
4. Within M.R the unbound metal binds to the E. coli with some rate constant α, and detaches
with some rate β (equilibrium process).
5. This process operates until the metal concentration is as low as desired.
3
2.2.3
Bioreactor Assumptions
• Our cell surface engineering technique has specificity to only one substance in the waste water
(i.e. the substance desired for removal).
• We will only observe exponential bacterial growth within the reactor, i.e. the biomass has been
cultured to sufficiently large optical density before introduction to bioreactor.
• Suspended metal particles attach to the surface of our E. coli at a rate proportional to the metal
concentration, and the ratio of unused area available on the surface.
• Bound metal particles detach from the surface of our E. coli back into the wastewater at a rate
proportional to the metal concentration bound to the E. coli.
• Once the E. coli that are pumped out of the main chamber end up in the flotation tank, they
are immediately removed - no more growth and no unbinding of metal occurs.
• We ignore all the spatial properties of pipes between the main reactor and flotation tank, as
soon as part of the solution leaves the main take, it immediately arrives in the flotation tank,
and vice versa.
• The pipe leading into the flotation tank sprinkles the solution over the flotation tank so lightly
that the E. coli cells do not sink into the tank - they float on the surface, allowing for immediate
removal.
• The solution within the main reactor is continuously mixed (can consider it to be a homogeneous
solution).
• All the E. coli in the flotation tank float.
• A higher dilution rate (D) does not affect the detachment rate.
2.3
Equations
The dynamics within the bioreactor are given by:
N
= Niµ(S) − N · (D + κ)
(1)
S
= −µ(S)N γ−1
(2)
)
(δN − Xw
Xu
= βXw − αXu ·
(3)
δN
)
(δN − Xw
Xw
= −βXw + αXu ·
− DXw
(4)
δN
Where, κ is the death rate, δ the maximum carrying capacity of the cells, α the binding rate of metal
to the cells and β the detachment rate of metal from the cells.
Equation (1) represents the change in E. coli cell concentration, equation (2) represents the change in
limiting substrate concentration, equation (3) represents the change in unbound metal concentration,
and equation (4) represents the change in concentration of metal bound to the surface of the E. coli.
In equation (3) we have βXw [9] which represents metal unbinding from the E. coli population,
it depends on β: the kinetic detachment rate, and Xw : the concentration of metal that is bound to
(δN−X
)
w
the surface of the E. coli. We are also removing metal from the solution by αXu ·
. This
δN
depends on the kinetic binding rate α, the concentration of unbound metal within the solution (Xu),
(δN−X
)
w
and the proportion of available binding sites on the E. coli, this is represented by
. Equation
δN
(4) is the negative of equation (3) but we also remove Xw at some dilution rate D.
4
2.4
Numerical Solution
Figure 2: Bioreactor Numerical Simulation
We solved the ODE system composed of equations (1) - (4) using the ODE45 solver in MATLAB (see
the end of the paper of the MATLAB code), using the initial conditions Ni = 0.4, S = 0.63, Xu =
0.5, Xw = 0 and parameters:
• γ = 1.1
• α = 1.5
• κ = 0.01
• µm = 0.8
• D=3
• Ks = 2.87 × 10−7
• β = 0.03
• δ = 1.5
The above simulation was done for removing Cobalt from mining wastewater effluent, it utilizes a
metal binding mechanism outlined in [2] which gave us a range of viable α. We got Ks and µm from
[8], β, and κ from [9]. We estimated reasonable values for D, γ. See §3.5 for calculating the value of
δ.
These results are incredibly promising, it shows that the bioreactor can theoretically be useful at
solving real problems, and that it operates within a reasonable amount of time (6.5 hours).
2.5
Sensitivity Analysis
In order to understand the dynamics within the bioreactor, we performed a sensitivity analysis that
estimates the relative effect of a single parameter on the performance of the system. Unfortunately
it is not possible to define an explicit function relating the different parameters to each other, we
5
must perform a ’naive’ analysis; we vary one parameter while holding all other constant, and see how
this change affects the performance of the system. Performance is measured by how long it takes the
bioreactor to remove the large majority of metal. We set some small ε = 10−4mol L−1 to define a low
threshold of acceptable metal concentration. For §2.5.1 - §2.5.6, we kept all the parameters at:
• γ = 1.1
• α = 1.5
• κ = 0.01
• µm = 0.8
• D=3
• Ks = 2.87 × 10−7
• β = 0.03
• δ = 1.5
and only varied the parameter who’s effect was measured.
2.5.1
The Effect of Dilution Rate on Bioreactor Performance
Figure 3: Dilution Rate (D) vs. Bioreactor Performance
(See Figure 3) The higher the dilution rate, the better the performance. Above dilution rate D = 1 it
slowly approaches an asymptote. We should not have a D that is too large because it is energy costly.
Also a larger D leads to faster deterioration of the bioreactor.
6
2.5.2
The Effect of Initial Metal Concentration on Bioreactor Performance
Figure 4: Initial Metal Concentration (Xu) vs. Bioreactor Performance
(See Figure 4) The higher the initial metal concentration within the bioreactor, the longer it takes to
operate.
2.5.3
The Effect of Attachment Rate on Bioreactor Performance
Figure 5: Attachment Rate (α) vs. Bioreactor Performance
7
(See Figure 5) The higher the attachment rate (α), the better the performance of the system. It does
approach an asymptote quite quickly though. This is because of the maximum carrying capacity, at
some point, there is just no more available space on the E. coli surface to bind more metal, no matter
the attachment rate.
2.5.4
The Effect of Detachment Rate on Bioreactor Performance
Figure 6: Detachment Rate (β) vs. Bioreactor Performance
(See Figure 6) The lower the detachment rate (β), the better the performance of the system.
2.5.5
The Effect of Initial Cellular Population on Bioreactor Performance
(See Figure 7) The higher the initial E. coli population (Ni), the better the system performs. After
some Ni ≈ 1, it approaches an asymptote.
8
Figure 7: Initial Cellular Population (Ni) vs. Bioreactor Performance
2.5.6
The Effect of Maximum Specific Growth Rate on Bioreactor Performance
Figure 8: Maximum Specific Growth Rate (µm) vs. Bioreactor Performance
(See Figure 8) The higher the maximum specific growth rate, the better the performance of the system.
2.6
Hypothetical Applications
The numerical simulations above were all based on parameters values specific to cobalt removal from
mining wastewater. If the E. coli are coupled with a binding mechanism that is not specific to cobalt,
9
then this generic bioreactor can be used for a large variety of applications (one simply has to change
α and β). One very promising application is the removal of pharmaceuticals, such as penicillin, from
municipal waste water. The reactor is viable for many different binding methods, strains of E. coli,
and initial concentrations of waste, as displayed in §2.5. This reactor is independent of volume, all
the parameters are relative to each other, this system would perform the same with 1 L of effluent as
with 1 × 1010 L of effluent (of course size of chambers would need to be scaled appropriately).
3
A Mechanistic Buoyancy Model: Maximum Carrying Ca-
pacity Determination
3.1
Goal
The buoyancy model has 2 main goals:
1. Determine the mean buoyant force experienced by the genetically engineered E. coli biomass
as a result of gas vesicle formation from ARG1 over-expression. This effectively determines a
mechanical upper bound for carrying capacity of the biomass.
2. Characterize and quantify flotation observations from wet lab experimentation. It is not well
known what the role of some of the secondary GVPs in the ARG1 construct is for gas vesicle
formation, and being able to quantify flotation facilities comparison of different combinations
of secondary GVPs for optimization of a gene construct specifically
engineered
for
biomass
flotation.
3.2
Mathematical Formulation
3.2.1
Concepts
• Newtonian mechanics:
F = ma ≡ mdv
t
-
F ≡ vector sum of forces on body
– m ≡ mass of body
- a ≡ acceleration as a function of time
- v ≡ velocity as a function of time
⃗
• Stokes-Einstein Drag:
F
D = −6πηRv
⃗
–
F
D ≡ drag force (note it opposes direction of velocity)
– η ≡ viscosity of fluid medium
- R ≡ (approximate) radius of body in spherical approximation.
• Integrating factor method to solve first order linear ODE:
dy
+ f(t)y = g(t)
(5)
dt
Let µ(t) = eβt. Multiply both sides by µ(t).
(
)
Chain rule
=⇒d
y(t)µ(t)
= h(t)µ(t)
dt
Integrate, divide by µ(t)
=⇒ y(t).
10
3.2.2
Assumptions
• The biomass separates into clumps that can be approximated by spheres.
• Buoyant force is constant.
• Cell motility and in-plane motion (motion perpendicular to the vertical axis) is negligible. This
effectively reduces our system to one dimension.
3.2.3
The Model
Stokes-Einstein, gravitational force near the surface of the earth Fg = mg (g ≈ 9.81 ms−2), Newton’s
second law lets us us write:
Fnet = Fg + FB + FD =⇒ mdv
= −mg + FB − 6πηRv
(6)
dt
Rearrange...
dv
6πηR
(FB
)
−
(7)
dt
m v=
m −g
(F
)
B
Simplify notation with α :=6πηR
g
we have:
m andβ:=
m −
dv
− αv = β
(8)
dt
Observe the model is now in the form of equation 5.
3.3
Analytic Solution
Now to get a solution we can directly apply the integrating factor method as described in the Concepts
sections.
The integrating factor is:
µ(t) = eαt
(9)
Multiplying both sides we have:
d
(
veαt) = βeαt
(10)
dt
Now we integrate, divide by integrating factor, and impose an initial condition of v(t = 0) = 0 ms−1
(not moving initially):
β(
(FB − mg)(
v(t) =
1−e−αt)≡
1−e−6mηRt)
(11)
α
6πηR
This is a closed form time-dependent solution for velocity.
Now we integrate one more time and impose an initial condition of z(t = 0) = 0 m:
)(
(FB − mg
m
(
z(t) =
t+
e6
m t −1))
(12)
gπηR
6πηR
Finally, we have a closed form time-dependent solution for vertical displacement.
NB: m, g, η, R are considered to be known, empirical constants. FB is left as a parameter that
can be determined in a curve fit regression to a time-series of vertical displacement data for a floating
biomass.
11
3.4
Determining Buoyant Force: Temporal Tracking Algorithm
Now that we have an expression for vertical displacement as a function of time which depends on
the buoyant force as a parameter. We want to determine the magnitude of this buoyant force. In
principle, this is not a difficult task as it can be solved using the built-in least-squares optimization
curve-fit functions, in say, MATLAB. The difficulty lies in acquiring experimental data to fit to. Our
proposed solution to this is to acquire images from a stationary point of view at evenly spaced, small
time intervals and perform image segmentation and analysis techniques to track the vertical position
of floating biomass frame-to-frame. In reality the biomass will be clumped into many clusters, so
clustering and labelling algorithms are employed. The advantage of using a visual tracking algorithm
to observe vertical displacement is that researchers can exactly characterize flotation, benchmark their
results, and calculate how different variables affect flotation.
(b) Binarized ROI
(a) Cultured E. coli expressing RFP
(c) Labelled ROI
Figure 9: Example of image segmentation and clustering
Below is pseudocode for a stochastic temporal tracking algorithm that is intended to maintain
consistent cluster labelling from frame-to-frame and account for clusters combining and splitting.
The time interval between frames should be chosen to be sufficiently small such that the probability
of two binding/un-binding events occurring between any two frames can be assumed to be zero.
• Load directory with images.
• For i up to number of frames (i.e. time steps):
- Read pixel data from image file.
- Manually set crop margins to ROI encapsulating the region of flotation.
(apply same
cropping margins in subsequent frames programatically).
- Convert from RGB to grayscale, perform thresholding using Otsu’s method.
12
- Binarize image based off of threshold, segmentation complete.
- Cluster and label binary image.
- Centroid clusters. Store labels j = 1, 2, 3, ..., N , positions (x(j),y(j)), and approximate
i
i
radius for each of the N objects identified in segmentation. Note the radius is updated at
each step to account for changes in cluster morphology affecting magnitude of drag force
at each time step.
- if i > 1:
∗ For all N labels j in frame i find position (x(j′ )
yi−1 s.t. the distance
i−1,yi−′ 1)withyi ≥
d((x(j′ )
i−1,yi−′ 1),(xij),yij)))isminimized.
∗ The label j in the ith frame is considered to be a child of the label j′ in the (i + 1)th
frame.
∗ End if-statement.
- End for-loop.
- Define function handle for solution to ODE model for buoyancy.
- Convert cluster tracking data from pixels to spacial units (based on pixel size of image),
curve-fit to each set of cluster branches.
- Take average of curve-parameter determination of FB for each branch, analyze distribution,
significance of fit, variance of mean.
But wait! We passively assumed that the number of clusters N′ in frame i − 1 was greater than the
number of clusters N in frame i. Not to worry, small modification...
1. N = N′:
• If this is the case, the logic holds and we simply have a child label j′ that is mapped to
parent label j.
2. N > N′:
• Physically, this corresponds to a cluster splitting.
• How do we handle this? No modification of the pseudo code algorithm is needed. There
will simply be two objects with labels j1, j2, respectively, to which some j′ is mapped to.
A branch in the tracking tree.
3. N < N′
• Non-trivial case.
• Switch search: For all N′ labels j′ in frame i + 1 find position (x(j),y(j)) with yi ≥ yi−1
i
i
s.t. the distance d((x(j′ )
i−1,yi−′ 1),(xij),yij)))isminimized.
• There will exist a label j in the ith frame such that two labels j′
1,j2 inthe(i+1)th frame
map to it. (2 clusters combine)
Note: The position tracking effectively forms a directed acyclic graph (DAG)
13
Figure 10: A few potential layers of a simple example with 3 initial clusters (nodes) to illustrate the
3 cases.
3.5
Discussion On Maximum Carrying Capacity
For this discussion, consider the buoyant force FB to be a known, constant value.
In order for the biomass to float the net force on it must be greater than zero. By equation 6
we have the following constraint on buoyant force:
FB ≥ mg + 6πηRv(t), ∀t
(13)
The right side of the above inequality has strictly positive monotonicity with respect to v. v is a
strictly increasing function so its max on a closed, bounded interval is at the upper bound. Let
vmax = v(tf ), where tf is the solution to z(tf ) = h, h being the height required from start point to
top of container, z being as given in equation 12. The inequality in equation 13 is thus equivalent to
the following:
FB ≥ mE.colig + 6πηRE.colivmax
(14)
Taking m to be the approximate mass of an E. coli cell (without bound particle), RE.coli to be the
approximate radius of an E. coli cell in a spherical approximation, this constraint on buoyant force
influences goals in the wet lab by setting a minimum requirement for observed buoyant force in mod-
ifying ARG1 for optimized flotation performance.
Once ARG1 is optimized for flotation and an estimate for FB is obtained the max carrying capacity
can be obtained as follows:
FB − 6πηRE.colivmax
m≤
(15)
g
Here m = mE.coli + mbound (sum of cell mass with mass of bound particle). Solving for mbound:
FB − 6πηRE.colivmax
mbound ≤
−mE.coli
(16)
g
14
Diving by the molecular or atomic mass, denote this mparticle of the substance to be sequestered we
have the following value for max carrying capacity δ in the bioreactor model:
1
[FB − 6πηRE.colivmax
]
δ=
−mE.coli
(17)
mparticle
g
But in the bioreactor model the effluent is pumped from the main reactor in such a way that cells will
not sink, so vmax can be set to 0 and we see:
1
[FB
]
δ≈
−mE.coli
(18)
mparticle g
See §2.3. This would influence the expression levels future cell surface engineering techniques. Too
many surface receptors could result in accumulation of too much particle mass, eliminating ability for
flotation.
4
Growth Dynamics
4.1
Objective
Our objective in this paper is to describe the growth dynamics of BL21 E. coli strain in Lysogeny
broth (LB). We have obtained two sets of observations from our wet lab team:
1.
25 observations with Δt = 30mins
2.
40 observations with Δt = 30mins
4.2
Approach
To model all the observations, we will use the Gompertz function[10], which can be expressed in the
following way:
f (t) = ae−be−ct
with:
b, c > 0
The reason for using this model is that it best describes the lag phase, the exponential phase and
stationary phase.
4.3
Results
The observations were given to us in OD600. We know that for bacterial cell cultures OD600 of
1.0 = 8.0 × 108cells/ml [6].
• For the first set of observations we fit the curve to our data.
15
8
10
8
Growth vs. time
Gompertz curve
7
6
5
4
3
2
1
0
0
5
10
15
20
25
time(30mins)
Figure 11: Fitted data
with the following parameters
General model:
ans(x)
= a∗exp(−b∗exp(−c∗x))
Coefficients
(with
95% c o n f i d e n c e bounds ) :
a =
3.876e+08
(2.782e+08,
4.97e+08)
b =
18.97
(−22.7,
60.64)
c
=
0.4938
(0.1035,
0.884)
• For the second set of observations we fit the same model. However, we obtain different param-
eters, since we enter the stationary phase after the the exponential phase, which we did not
observe in the first data set. But we also notice that they fall into the 95% confidence interval
for the parameters of the first observations. Hence we can assume, with some certainty, that the
parameters we got for this curve fitting process would also fit for the previous one.
16
8
10
7
cells/ml vs. time
Fitted Curve
6
5
4
3
2
1
0
0
5
10
15
20
25
time(30mins)
Figure 12: Fitted data
The parameters are as follows:
General model:
ans(x)
= a∗exp(−b∗exp(−c∗x))
Coefficients
(with
95% c o n f i d e n c e bounds ) :
a =
5.707e+08
(5.355e+08,
6.059e+08)
b =
5.127
(3.785,
6.469)
c
=
0.1989
(0.1603,
0.2376)
5
Discussion and Future Directions
Based on the results of our differential bioreactor model we postulate that a bioreactor of this design
could perform at appropriately small time-scales with a sufficiently optimized flotation construct in a
bioremediation context. This bodes well for future laboratory endeavours where the bioreactor schema
along with an engineered cell-line optimized for flotation from gas vesicle formation could be tested
in a small scale laboratory model of the system to test its empirical performance. This is useful for
model validation, and proof of concept.
An improvement to the bioreactor could be seen in the form of incorporation of an experimental
growth dynamics model as demonstrated in the Growth Dynamics section with the Goempertz equa-
tion. The Monod equation was chosen initially for this model because it is scalable and its parameters
are known for BL21; however, it only describes the exponential phase of growth and it is not guar-
anteed that the biomass will remain in exponential phase during the bioreactor process. Another
shortcoming is that it is a single-substrate model and in practice this would not be the case as the
likely nutrient source would be from carbohydrates, lipids, and organics present in the waste water
17
effluent. For future application-based analysis using the bioreactor model, Goempertz coefficients
could be determined for a biomass of industrially relevant size for the volume demands dictated by
the industry requiring bioremediation using a similar experimentation and analysis technique as de-
scribed in the Growth Dynamics section.
To further enrich the analysis of the behaviour of the bioreactor model, in addition to the exist-
ing sensitivity analysis, identification and characterization of stable points resulting from different
combinations of parameter values could be performed.
Application of the stochastic temporal tracking algorithm to real flotation data (images) to estimate
the buoyant force for ARG1 and compare different modifications of ARG1 to determine an optimal
gene combination for flotation.
A next major step is coupling cellular flotation assays with cell surface engineering for application-
based testing for different bioremediation tasks. Much research has been done relating to selecting
peptide sequences for metal binding [4][5]. Existing techniques could be coupled with a signal trans-
duction pathway to induce expression of ARG1 or modified versions of it optimized for flotation.
18
6
Appendix
Nomenclature
α
Attachment Rate (h−1)
β
Detachment Rate (h−1)
δ
Maximum Carrying Capacity Constant (h−1)
η
Dynamic viscosity of waste water (kg ms−1)
γ
Yield Constant (unitless)
κ
E. coli Death Rate (h−1)
µm Maximum Specific Growth Rate (h−1)
D=F
V DilutionRate(h−1)
F
Flow Rate (Lh−1)
Ks Monod Constant (mol L−1)
mE.coli Aproximate mass of E. coli cell
mparticle Atomic or molecular mass of substance to be sequestered (g mol−1)
Ni
Initial Cellular Concentration (mol L−1)
RE.coli Approximate radius of E. coli cell in spherical approximation (m)
Si
Initial Limiting Substrate Concentration (mol L
−1)
V
Total Volume (L)
vmax Maximum achievable velocity based on vertical displacement specification.
Xui Initial Metal Concentration within Waste water (mol L−1)
Xw Concentration of Bound Metal (mol L−1)
7
References
[1] Oladeji, S. O., & Saeed, M. D. (2015). Assessment of cobalt levels in wastewater, soil and vegetable
samples grown along Kubanni stream channels in Zaria, Kaduna State, Nigeria. African Journal
of Environmental Science and Technology, 9(10), 765-772. doi:10.5897/ajest2015.1969
[2] Whittaker MM, Mizuno K, Bächinger HP, Whittaker JW. Kinetic Analysis of the Metal
Binding Mechanism of Escherichia coli Manganese Superoxide Dismutase. Biophysical Journal.
2006;90(2):598-607. doi:10.1529/biophysj.105.071308.
[3] Monod substrate affinity constant
(Ks) for strain ML30
of E. coli growing on glu-
cose as the only source of carbon and energy.
(n.d.). Retrieved October
13,
2018
from
[4] N., Cetinel, S., Omar, S. I., Tuszynski, J. A., & Montemagno, C. (2017). A computational method
for selecting short peptide sequences for inorganic material binding. Proteins: Structure, Function,
and Bioinformatics, 85(11), 2024-2035. doi:10.1002/prot.25356
19
[5] Li, Pengsong & Tao, Huchun. (2013). Cell surface engineering of microorganisms towards adsorp-
tion of heavy metals. Critical reviews in microbiology. 41. 10.3109/1040841X.2013.813898.
[6] Bacterial
cell
number
(OD600)
(n.d.).
Retrieved October
13,
2018
from
[7] Bengt Carlsson.
(2009). An introduction to modeling of bioreactors. Dept of Systems and
Control, Information Technology Uppsala University. Retrieved October
13,
2018
from
[8] Monod substrate affinity constant (Ks) for strain ML30 of E. coli growing on glucose as the only
source of carbon and energy. (n.d.). Retrieved October 13, 2018
[9] Chen, L., & Chai, L.
(2005). Mathematical Model and Mechanisms for Biofilm Wastewa-
ter Treatment Systems. World Journal of Microbiology and Biotechnology, 21(8-9), 1455-1460.
doi:10.1007/s11274-005-6565-2
[10] Parolini, N, & Carcano, S. (2009). A MODEL FOR CELL GROWTH IN BATCH BIOREAC-
TORS. Faculty of Systems Engineering, Polytechnic University of Milan. Retrieved October 13,
2018 from https://www.politesi.polimi.it/bitstream/10589/2082/1/2010 07 Carcano.pdf.
20
The MATLAB Bioreactor ODE Code
gamma = 1 . 1 ; % y i e l d c o n s t a n t
(?)
kappa = 0.01;
% death rate hˆ−1
D = 3; % dilution /flow rate hˆ−1
beta = 0.03;
% detachment rate of particle from biomass surface hˆ−1
alpha
= 1.5;
% attachment rate hˆ−1
mu m = 0 . 8 ; % maximum s p e c i f i c g r o w t h r a t e
(monod equation ) u n i t l e s s
K s = 2.87∗10ˆ( −7); % monod constont mol Lˆ−1
delta
= 1.5;
% maximum c a r r y i n g c a p a c i t y
N i = 0.4; % i n i t i a l cellular concentration Lˆ−1
S i
= 0.63;
% i n i t i a l substrate concentration Lˆ−1
Xu i = 0.5;
% i n i t i a l unbound metal concentration Lˆ−1
Xw i = 0; % i n i t i a l bound metal concentration Lˆ−1
epsilon
= 0.0001; % Acceptable metal concentration threshold
f
= @( t , x)
[N i
∗
((mu m ∗ x(2))
/
(K s + x(2)))
− x(1)
∗
(D + kappa ) ;
−((mu m ∗ x (2))
/
(K s + x(2)))
∗ x(1)
/ gamma;
beta
∗ x(4)
− alpha
∗ x(3)
∗
(delta
∗ x(1)
− x(4))/(delta
∗ x(1));
−(beta ∗ x(4) − alpha
∗ x(3)
∗
(delta
∗ x(1)
− x(4))/(delta
∗ x(1)))
− D∗x ( 4 ) ] ;
[t,xa]
= ode45( f , [ 0
50],[N i S i Xu i Xw i]);
fprintf(’ode solved\n’)
ind
= find(xa (: ,3) <= epsilon );
xap = xa (1: ind , : ) ;
tp
= t(1:ind);
figure;
plot(tp,xap(: ,1));
% Cell concetration mol Lˆ−1
hold on
plot(tp,xap(: ,2));
% Substrate concentration
hold on
plot(tp,xap(: ,3))
% Unbound metal concentration mol Lˆ−1
hold on
plot(tp,xap(: ,4))
% Bound metal concentration mol Lˆ−1
hold on
set(findall(gca,
’Type ’ ,
’Line’), ’LineWidth’ ,2);
legend( ’E. c o l i Pop ’ ,
’Substrate ’,’Unbound Particle ’,
’Bound p a r t i c l e in
MR’ )
xlabel(’Time
(h)’), ylabel(’mol/L’)
21