September 25, 2017, Monday
University of Colorado at Boulder Search A to Z Campus Map University of Colorado at Boulder CU 
Search Links


MBW:Flocculation Dynamics - a PDE Model

From MathBio

Jump to: navigation, search

Executive Summary

In "Klebsiella pneumoniae Flocculation Dynamics", Bortz, et. al [1] propose a nonlinear transport partial differential equation (PDE) as a model for the flocculation dynamics of bacterial aggregates in suspension. One goal for this Wiki page is to describe their model. Another goal is to provide enhanced detail on how the authors approached solving the PDE numerically. Finally, the most important goal is to present and explain a Matlab program that approximates a solution to the PDE.

Overview

  1. Mathematics Used - The following model heavily relies on Partial Differential Equations to model aggregation, flocculation, and growth (respectively) within a closed system (i.e. a body). Using the Method of Lines (MOL) technique, the governing PDE is reduced -through discretization of the spatial dimensions- into a system on Ordinary Differential Equations. The system of ODE's can then be solved in numerous ways. For other examples of PDEs in mathematical biology, look at Traveling Waves in Excitable Media and Gravitational Effects on Blood Flow.
  2. Type of Model - This model attempts to obtain numerical rates for aggregation, fragmentation and growth of clusters of bacteria.
  3. Biological System Studied- Much research is being done to study the interaction of miscro-organisms, inside the body and outside. This research can have applications from epidemiology to water treatment to evolutionary biology and the origination of multicellular life.

Context/Biological Phenomenon under Consideration

"Klebsiella pneumoniae" are a common source of infections in the bloodstream. How these bacteria behave in the bloodstream, and in particular, how they behave in suspension (as opposed to how they behave attached to a substrate), provides a focus for the paper mentioned above and for the study in this Wiki. We consider three factors in the flocculation dynamics of the bacteria, aggregation, fragmentation, and growth. Aggregation simply means an individual bacteria or group of bacteria joining together with another individual bacteria or group of bacteria to form one larger group (or aggregate). Fragmentation is the reverse process, where an aggregate of bacteria breaks into smaller aggregates, meaning each new aggregate consists of fewer bacteria than the original clump. Finally, growth refers to the increase in size (volume) of an aggregate.

  • For those of you food lovers out there: Believe it or not, flocculation has a large impact in Cheese Production

Mathematical Model with Variable/Parameter Definitions

Bortz, et al., propose the following governing PDE as a model, and solving it provides a focus for this wiki

b_{t}+(G(x)b)_{x}=A(x,b)+F(x,b)

with G(\underline {x})b(t,\underline {x})=0 and b(0,x)=b_{0}(x). Figure 1 shows a representation of the processes included in this model.

Figure 1: Aggregation, Fragmentation, and Growth Schematic.
Variables
  • x - volume of bacterial aggregate (in fL)
  • \underline {x}=2fL, \overline {x}=1000fL - assumed minimum and maximum volumes
  • t - time (in minutes unless otherwise stated)
  • b(t,x) - density function representing the number of aggregates of size x per fL
  • G(x)=\gamma _{G}x(1-{\frac  {x}{\overline {x}}}) - assumed logistic growth rate of bacterial aggregates of size x
  • A(x,b)=A_{{in}}-A_{{out}} - net rate at which bacterial flocs aggregate
    • A_{{in}}={\frac  {1}{2}}\int _{{\underline {x}}}^{{x-\underline {x}}}K_{A}(y,x-y)b(t,y)b(t,x-y)dy for x\in [2\underline {x},\overline {x}]
    • A_{{out}}=b(t,x)\int _{{\underline {x}}}^{{\overline {x}-x}}K_{A}(x,y)b(t,y)dy for x\in [\underline {x},\overline {x}-\underline {x}]
    • K_{A}(x,y)=\gamma _{A}\left({\frac  {\epsilon }{\nu }}\right)^{{1/2}}\left(x^{{1/3}}+y^{{1/3}}\right)^{3} - assumed turbulent mixing kernel describing the rate at which flocs of volumes x and y join to make a floc of size x + y
  • F(x,b)=F_{{in}}-F_{{out}} - net rate at which bacterial flocs fragment
    • F_{{in}}=\int _{x}^{{\overline {x}}}\Gamma (x,y)K_{F}(y)b(t,y)dy for x\in [\underline {x},\overline {x}-\underline {x}]
    • F_{{out}}={\frac  {1}{2}}K_{F}(x)b(t,x)\triangle x for x\in [2\underline {x},\overline {x}]
    • K_{F}(x)=\gamma _{F}(x-\underline {x})^{{1/3}} - describes rate at which flocs fragment
    • \Gamma (x,y)={\begin{cases}3/y;&{\frac  {y}{3}}<x\leq {\frac  {2y}{3}},\\0;&{\textrm  {otherwise.}}\end{cases}} - describes probability density of daughter flocs for the fragmentation of a parent floc of size y
Parameters
  • \gamma _{G}=6.8\times 10^{{-4}}/{\textrm  {min}}
  • \gamma _{A}=2.7\times 10^{{-15}}/{\textrm  {fL}}^{2}
  • \gamma _{F}=6.6\times 10^{{-5}}/(\mu {\textrm  {mfLmin)}}
  • \epsilon =4.43\times 10^{{-4}}{\textrm  {m}}^{2}/{\textrm  {s}}^{3}
  • \nu =1.99\times 10^{{-6}}{\textrm  {m}}^{2}/{\textrm  {s}}

Analysis (Analytical and Computational)

  • In order to solve the governing PDE, we will employ a strategy known as Method of Lines (MOL). In general, to execute the MOL technique, we discretize our PDE in its spatial variable(s), which will result in a system of ODEs. We can then use our favorite ODE solver to advance the PDE solution in discrete time steps.

Discretization Scheme

Bortz, et al., provide a detailed analysis of the discretization scheme and a proof for the convergence of the scheme including the allowable space of functions, H, we can use for our density function. An assumption we make for the rest of this Wiki is that of an equispaced mesh. To produce the scheme, we need a set of basis elements for i=1,\dots ,N,

\beta _{i}^{N}(x)={\begin{cases}1;\;x_{{i-1}}^{N}\leq x<x_{i}^{N}\\0;\;{\textrm  {otherwise}}\end{cases}}

where x_{i} represents a node in the mesh from the minimum volume x_{0} to the maximum volume x_{N}. These functions form an orthogonal basis for

H^{N}=\left\{h\in H|\;h=\sum _{{i=1}}^{N}\alpha _{i}\beta _{i}^{N},\;\alpha \in {\mathbb  {R}}\right\}

with projections \pi ^{N}:H\rightarrow H^{N}

\pi ^{N}h=\sum _{{j=1}}^{N}\alpha _{j}\beta _{j}^{N};{\textrm  {where}}\alpha _{j}={\frac  {1}{\triangle x}}\int _{{x_{{j-1}}^{N}}}^{{x_{j}^{N}}}h(x)dx.

Finally, we need one more operator, {\mathcal  {G}}^{N}, defined as

({\mathcal  {G}}^{N}h)(x)=\sum _{{j=1}}^{N}{\frac  {1}{\triangle x}}\left(G(x_{{j-1}})\alpha _{{j-1}}-G(x_{{j}})\alpha _{{j}}\right)\beta _{j}^{N}(x).

All of this provides the spatial discretization needed for MOL to provide the resulting system of ODEs

b_{t}^{N}={\mathcal  {G}}^{N}b^{N}+\pi ^{N}({\mathcal  {A}}(b^{N})+{\mathcal  {F}}(b^{N}))

b^{N}(0,x)=\pi ^{N}b_{0}(x).

Each of the operators in the discretized system generate a matrix (or vector) that in combination provide the framework for an ODE solver. The following gives a detailed description of each term in the discretized system. We also include the code associated with each.

Initial Conditions

Initial Conditions and the left hand side vector, [b^{N}]_{t}, in our ODE system.

Based on actual data taken five minutes after the biofilm colony was placed in the swirling flask, the authors determined a bi-exponential initial density function. The function used in this Wiki has been updated since the publishing of the paper.

b_{0}(x)=89.828e^{{-0.0035303x}}+1946.6\times 10^{6}e^{{-1.3159x}}

The projection into H^{N} of these initial conditions is derived as follows: \pi ^{N}(b_{0}(x))=\sum _{{j=1}}^{N}\alpha _{j}\beta _{j}^{N}(x){\textrm  {where}}\alpha _{j}={\frac  {1}{\triangle x}}\int _{{x_{{j-1}}^{N}}}^{{x_{j}^{N}}}b_{0}(x)dx.

So we have, \alpha _{1}={\frac  {1}{\triangle x}}\int _{{x_{0}^{N}}}^{{x_{1}^{N}}}b_{0}(x)dx,\;\alpha _{2}={\frac  {1}{\triangle x}}\int _{{x_{1}^{N}}}^{{x_{2}^{N}}}b_{0}(x)dx,{\textrm  {etc.}}

Then to develop the left hand side of our system of differential equations recall that \beta _{j}^{N}=1 when x_{{j-1}}\leq x<x_{j}. So for instance, \beta _{1}(x_{0})=1, but \beta _{1}(x_{i})=0 for any i\neq 1. Now we have

\pi ^{N}b_{0}(x_{0})=\alpha _{1}\beta _{1}^{N}(x_{0})+\alpha _{2}\beta _{2}^{N}(x_{0})+\dots +\alpha _{N}\beta _{N}^{N}(x_{0})=\alpha _{1}

\pi ^{N}b_{0}(x_{1})=\alpha _{1}\beta _{1}^{N}(x_{1})+\alpha _{2}\beta _{2}^{N}(x_{1})+\dots +\alpha _{N}\beta _{N}^{N}(x_{1})=\alpha _{2}

\vdots

\pi ^{N}b_{0}(x_{N})=\alpha _{1}\beta _{1}^{N}(x_{N})+\alpha _{2}\beta _{2}^{N}(x_{N})+\dots +\alpha _{N}\beta _{N}^{N}(x_{N})=\alpha _{N}

In other words,

[b^{N}]_{t}={\begin{bmatrix}\alpha _{1}\\\alpha _{2}\\\vdots \\\alpha _{N}\end{bmatrix}}_{t}

To implement this in Matlab and actually form our initial conditions on the mesh, we wrote the following function. Of note, the first iteration of this function used symbolic integration on the original b_{0} defined above for each calculation on the mesh. On a suggestion from Dr. Bortz, we ran a symbolic integration on b_{0} using symbolic variables x_{1} and x_{2} to generate the expression for the code below. With a mesh of 50 increments from x_{0} to x_{N}, this technique sped up our code significantly.

function c = intinitcond2(x1,x2) % This function will integrate the symbolic initial function from volume % x1 to volume x2.

c = 19466000000000/(13159*exp((13159*x1)/10000)) - ...

   19466000000000/(13159*exp((13159*x2)/10000)) + ...
   809098694654873829376/(31798115529012125*...
   exp((254384924232097*x1)/72057594037927936)) - ...
   809098694654873829376/(31798115529012125*...
   exp((254384924232097*x2)/72057594037927936)) ;

Additionally, we can see the general pattern for the left hand side vector, b^{N}]_{t}, in our ODE system.

{\begin{bmatrix}\sum _{{i=1}}^{N}\alpha _{i}\beta _{i}^{N}(x_{0})\\\sum _{{i=1}}^{N}\alpha _{i}\beta _{i}^{N}(x_{1})\\\vdots \\\sum _{{i=1}}^{N}\alpha _{i}\beta _{i}^{N}(x_{{N-1}})\end{bmatrix}}_{t}

Approximating generators

{\mathcal  {G}}^{N}:H^{N}\rightarrow H^{N} of the infinitesimal generator

This was defined as

({\mathcal  {G}}^{N}h)(x)=\sum _{{j=1}}^{N}{\frac  {1}{\triangle x}}\left(G(x_{{j-1}})\alpha _{{j-1}}-G(x_{{j}})\alpha _{{j}}\right)\beta _{J}^{N}(x),

for h\in H^{N}. Recall, G(x) represents an assumed logistical growth with maximum size of 1000 fL, where G(x)=\gamma _{G}x(1-{\frac  {x}{\overline {x}}}). When the approximating growth generator is applied to the projected density function, b^{N}, in our system of ODEs, we get the relation

{\mathcal  {G}}^{N}\beta _{J}^{N}={\frac  {1}{\triangle x}}G(x_{j}^{N})\beta _{{j+1}}^{N}-{\frac  {1}{\triangle x}}G(x_{j}^{N})\beta _{j}^{N}.

To see why this is true, consider application of the approximating generator to \beta _{2}^{N}, and note the substitution, \alpha _{j}={\frac  {1}{\triangle x}}\int _{{x_{{j-1}}^{N}}}^{{x_{j}^{N}}}\beta _{2}^{N}(x)dx. A key point to remember below is that \beta _{2}^{N}(x)=1 from x_{1} to x_{2} and it is zero everywhere else. So when we integrate it across those nodes we get \triangle x, whereas it integrates to zero across any other pair of nodes.

{\mathcal  {G}}^{N}\beta _{2}^{N}(x)=\sum _{{j=1}}^{N}{\frac  {1}{\triangle x}}\left(G(x_{{j-1}}){\frac  {1}{\triangle x}}\int _{{x_{{j-2}}^{N}}}^{{x_{{j-1}}^{N}}}\beta _{2}^{N}(x)dx-G(x_{j}){\frac  {1}{\triangle x}}\int _{{x_{{j-1}}^{N}}}^{{x_{j}^{N}}}\beta _{2}^{N}(x)dx\right)\beta _{J}^{N}(x)

={\frac  {1}{\triangle x}}\left\{\left[G(x_{0}){\frac  {1}{\triangle x}}\int _{{x_{{-1}}^{N}}}^{{x_{0}^{N}}}\beta _{2}^{N}(x)dx-G(x_{1}){\frac  {1}{\triangle x}}\int _{{x_{0}^{N}}}^{{x_{1}^{N}}}\beta _{2}^{N}(x)dx\right]\beta _{1}^{N}(x)\right.

+\left[G(x_{1}){\frac  {1}{\triangle x}}\int _{{x_{0}^{N}}}^{{x_{1}^{N}}}\beta _{2}^{N}(x)dx-G(x_{2}){\frac  {1}{\triangle x}}\int _{{x_{1}^{N}}}^{{x_{2}^{N}}}\beta _{2}^{N}(x)dx\right]\beta _{2}^{N}(x)

\left.+\left[G(x_{2}){\frac  {1}{\triangle x}}\int _{{x_{1}^{N}}}^{{x_{2}^{N}}}\beta _{2}^{N}(x)dx-G(x_{3}){\frac  {1}{\triangle x}}\int _{{x_{2}^{N}}}^{{x_{3}^{N}}}\beta _{2}^{N}(x)dx\right]\beta _{3}^{N}(x)+\cdots \right\}

={\frac  {1}{\triangle x}}\left\{0-0+0-G(x_{2}){\frac  {\triangle x}{\triangle x}}\beta _{2}^{N}(x)+G(x_{2}){\frac  {\triangle x}{\triangle x}}\beta _{3}^{N}(x)-0+0\cdots \right\}

={\frac  {1}{\triangle x}}G(x_{2})\beta _{3}^{N}(x)-{\frac  {1}{\triangle x}}G(x_{2})\beta _{2}^{N}(x)

Once we have the general relationship above, we can determine the contributions to each component of the left hand side of our ODE system. To illustrate, consider

{\mathcal  {G}}^{N}\beta _{1}^{N}(x)={\frac  {1}{\triangle x}}G(x_{1})\beta _{2}^{N}(x)-{\frac  {1}{\triangle x}}G(x_{1})\beta _{1}^{N}(x)

therefore {\mathcal  {G}}^{N}\beta _{1}^{N}(x_{0})={\frac  {1}{\triangle x}}G(x_{1})\beta _{2}^{N}(x_{0})-{\frac  {1}{\triangle x}}G(x_{1})\beta _{1}^{N}(x_{0})=0-{\frac  {1}{\triangle x}}G(x_{1})

{\mathcal  {G}}^{N}\beta _{1}^{N}(x_{1})={\frac  {1}{\triangle x}}G(x_{1})\beta _{2}^{N}(x_{1})-{\frac  {1}{\triangle x}}G(x_{1})\beta _{1}^{N}(x_{1})={\frac  {1}{\triangle x}}G(x_{1})-0

{\mathcal  {G}}^{N}\beta _{1}^{N}(x_{2})={\frac  {1}{\triangle x}}G(x_{1})\beta _{2}^{N}(x_{2})-{\frac  {1}{\triangle x}}G(x_{1})\beta _{1}^{N}(x_{2})=0={\mathcal  {G}}^{N}\beta _{1}^{N}(x_{3},x_{4},\dots x_{N})

This shows that in row one of our vector for $\mathcal{G}^N b^N$ we have the negative contribution from \beta _{1}^{N} of {\frac  {1}{\triangle x}}G(x_{1}), and in row two, \beta _{1}^{N} contributes {\frac  {1}{\triangle x}}G(x_{1}). Let's consider the contributions from \beta _{2}^{N} which should clearly illustrate the pattern for the rest of the vector for {\mathcal  {G}}^{N}b^{N}.

{\mathcal  {G}}^{N}\beta _{2}^{N}(x)={\frac  {1}{\triangle x}}G(x_{2})\beta _{3}^{N}(x)-{\frac  {1}{\triangle x}}G(x_{2})\beta _{2}^{N}(x)

therefore

{\mathcal  {G}}^{N}\beta _{2}^{N}(x_{0})={\frac  {1}{\triangle x}}G(x_{2})\beta _{3}^{N}(x_{0})-{\frac  {1}{\triangle x}}G(x_{2})\beta _{2}^{N}(x_{0})=0

{\mathcal  {G}}^{N}\beta _{2}^{N}(x_{1})={\frac  {1}{\triangle x}}G(x_{2})\beta _{3}^{N}(x_{1})-{\frac  {1}{\triangle x}}G(x_{2})\beta _{2}^{N}(x_{1})=0-{\frac  {1}{\triangle x}}G(x_{2})

{\mathcal  {G}}^{N}\beta _{2}^{N}(x_{2})={\frac  {1}{\triangle x}}G(x_{2})\beta _{3}^{N}(x_{2})-{\frac  {1}{\triangle x}}G(x_{2})\beta _{2}^{N}(x_{2})={\frac  {1}{\triangle x}}G(x_{2})-0

{\mathcal  {G}}^{N}\beta _{2}^{N}(x_{3})={\frac  {1}{\triangle x}}G(x_{2})\beta _{3}^{N}(x_{3})-{\frac  {1}{\triangle x}}G(x_{2})\beta _{2}^{N}(x_{2})=0={\mathcal  {G}}^{N}\beta _{2}^{N}(x_{4},x_{5},\dots x_{N})

The pattern then emerges as follows,

[{\mathcal  {G}}^{N}]b^{N}={\begin{bmatrix}-{\frac  {1}{\triangle x}}G(x_{1})\ &0&\cdots &\cdots &0\\{\frac  {1}{\triangle x}}G(x_{1})&-{\frac  {1}{\triangle x}}G(x_{2})&0&\cdots &0\\\vdots \\0&\cdots &0&{\frac  {1}{\triangle x}}G(x_{{N-1}})&-{\frac  {1}{\triangle x}}G(x_{N})\end{bmatrix}}b^{N}

The code to build this matrix follows.

% Build Growth Matrix

gammaG = 6.8e-4; % - parameter given

G = @(xx) gammaG.*xx.*(1-(xx./xmax));

Growthvec = zeros(N,N);

Growthvec(1,1) = -G(x(2));

for r = 2:N

   Growthvec(r,r-1) = G(x(r));
   Growthvec(r,r) = - G(x(r+1));

end

% The whole matrix has a factor of 1/delx

Growthvec = Growthvec./delx;

Fragmentation

Projected fragmentation matrix, \pi ^{N}({\mathcal  {F}}(b^{N}))

K_{F}(x) represents the fragmentation kernel. For this analysis, we used the kernel described by the authors as K_{F}(x)=\gamma _{F}(x-\underline {x})^{{1/3}}. Furthermore \Gamma (x,y) represents the authors' assumed uniform post-fragmetation distribution, centered around y/2, with width, y/3

\Gamma (x,y)={\begin{cases}3/y;&{\frac  {3}{y}}<x\leq {\frac  {2y}{3}},\\0;&{\textrm  {otherwise.}}\end{cases}}

To generate the probability distribution we built the following MATLAB function

function g = probdist(x1,x2)

% Determines the value of the post frag distribution function

% for inputs, x1 and x2, where x1 assumed less than x2

if (x1 > x2/3) && (x1 <= 2*x2/3)

   g = 3/x1;

else

   g =0;

end

We can then generate \pi ^{N}{\mathcal  {F}}(b^{N})

[\pi ^{N}{\mathcal  {F}}(b^{N})]={\begin{bmatrix}\sum _{{j=2}}^{N}\Gamma (x_{0},x_{{j-1}})K_{F}(x_{{j-1}})\alpha _{j}\triangle x\\-{\frac  {1}{2}}K_{F}(x_{1})\alpha _{2}\triangle x+\sum _{{j=3}}^{N}\Gamma (x_{1},x_{{j-1}})K_{F}(x_{{j-1}})\alpha _{j}\triangle x\\\vdots \\-{\frac  {1}{2}}K_{F}(x_{{N-2}})\alpha _{{N-1}}\triangle x+\Gamma (x_{{N-2}},x_{{N-1}})K_{F}(x_{{N-1}})\alpha _{N}\triangle x\\-{\frac  {1}{2}}K_{F}(x_{{N-1}})\alpha _{N}\triangle x\end{bmatrix}}

To see how this matrix develops from the actual model, we will demonstrate an example where we set N=4. Our partition then starts at x_{0}=2fL and ends at x_{4}=1000fL. The fragmentation operator acting on b^{4} is

{\mathcal  {F}}(b^{4})(x)=\int _{x}^{{x_{4}}}\Gamma (x,y)K_{F}(y)b^{4}(y)dy{\textrm  {;}}\;x\in [x_{0},x_{4}-x_{0}]

-{\frac  {1}{2}}K_{F}(x)b^{4}(x)\triangle x{\textrm  {;}}\;x\in [2x_{0},x_{4}]

and its projection is defined as

\pi ^{4}({\mathcal  {F}}(b^{4}))=\sum _{{j=1}}^{4}\gamma _{j}\beta _{j}^{4}(x){\textrm  {where}}\gamma _{j}=\int _{{x_{{j-1}}}}^{{x_{j}}}{\mathcal  {F}}(b^{4})(x)dx

Now consider row one of [\pi ^{N}{\mathcal  {F}}(b^{N})], so

\pi ^{N}{\mathcal  {F}}(b^{N})(x_{0})=\sum _{{j=1}}^{4}\gamma _{j}\beta _{j}^{4}(x_{0})=\gamma _{1}

where,

\gamma _{1}={\frac  {1}{\triangle x}}\int _{{x_{0}}}^{{x_{1}}}\left[\int _{{x_{0}}}^{{x_{4}}}\Gamma (x_{0},y)K_{F}(y)b^{4}(y)dy{\textrm  {;}}\;x\in [x_{0},x_{4}-x_{0}]\right]dx-0{\textrm  {since}}x_{0}\notin [2x_{0},x_{4}]

={\frac  {1}{\triangle x}}\int _{{x_{0}}}^{{x_{1}}}\left[\int _{{x_{0}}}^{{x_{1}}}\Gamma (x_{0},y)K_{F}(y)b^{4}(y)dy+\int _{{x_{1}}}^{{x_{2}}}\Gamma (x_{0},y)K_{F}(y)b^{4}(y)dy\right.

\left.+\int _{{x_{2}}}^{{x_{3}}}\Gamma (x_{0},y)K_{F}(y)b^{4}(y)dy+\int _{{x_{3}}}^{{x_{4}}}\Gamma (x_{0},y)K_{F}(y)b^{4}(y)dy\right]dx

For each of the inter integrals we'll use a right hand estimate. Also, note that b^{4}(x)=\sum _{{j=1}}^{4}\alpha _{j}\beta _{J}^{4}(x) and our basis elements \beta _{J}^{4}(x)=1 for the prior node, i.e., \beta _{1}^{4}(x_{0})=1 otherwise it is zero. Then

\gamma _{1}={\frac  {1}{\triangle x}}\int _{{x_{0}}}^{{x_{1}}}\left[\Gamma (x_{0},x_{1})K_{F}(x_{1})\alpha _{2}\triangle x+\Gamma (x_{0},x_{2})K_{F}(x_{2})\alpha _{3}\triangle x\right.

\left.+\Gamma (x_{0},x_{3})K_{F}(x_{3})\alpha _{4}\triangle x+\Gamma (x_{0},x_{4})K_{F}(x_{4})\times 0\times \triangle x\right]dx

=\Gamma (x_{0},x_{1})K_{F}(x_{1})\alpha _{2}\triangle x+\Gamma (x_{0},x_{2})K_{F}(x_{2})\alpha _{3}\triangle x+\Gamma (x_{0},x_{3})K_{F}(x_{3})\alpha _{4}\triangle x

=\sum _{{j=2}}^{4}\Gamma (x_{0},x_{{j-1}})K_{F}(x_{{j-1}})\alpha _{j}\triangle x

Now consider row two where \pi ^{N}{\mathcal  {F}}(b^{N})(x_{1})=\gamma _{2}, and

\gamma _{2}={\frac  {1}{\triangle x}}\int _{{x_{1}}}^{{x_{2}}}\left[\int _{{x_{1}}}^{{x_{4}}}\Gamma (x_{1},y)K_{F}(y)b^{4}(y)dy\right]dx-{\frac  {1}{\triangle x}}\int _{{x_{1}}}^{{x_{2}}}{\frac  {1}{2}}K_{F}(x_{1})b^{4}(x_{1})\triangle xdx

={\frac  {1}{\triangle x}}\int _{{x_{1}}}^{{x_{2}}}\left[\int _{{x_{1}}}^{{x_{2}}}\Gamma (x_{1},y)K_{F}(y)b^{4}(y)dy+\int _{{x_{2}}}^{{x_{3}}}\Gamma (x_{1},y)K_{F}(y)b^{4}(y)dy\right.

\left.+\int _{{x_{3}}}^{{x_{4}}}\Gamma (x_{1},y)K_{F}(y)b^{4}(y)dy\right]dx-{\frac  {1}{2}}K_{F}(x_{1})\alpha _{2}\triangle x

=-{\frac  {1}{2}}K_{F}(x_{1})\alpha _{2}\triangle x+{\frac  {1}{\triangle x}}\int _{{x_{1}}}^{{x_{2}}}\left[\Gamma (x_{1},x_{2})K_{F}(x_{2})b^{4}(x_{2})\triangle x+\Gamma (x_{1},x_{3})K_{F}(x_{3})b^{4}(x_{3})\triangle x\right.

\left.+\Gamma (x_{1},x_{4})K_{F}(x_{4})b^{4}(x_{4})\triangle x\right]dx

=-{\frac  {1}{2}}K_{F}(x_{1})\alpha _{2}\triangle x+\Gamma (x_{1},x_{2})K_{F}(x_{2})\alpha _{3}\triangle x+\Gamma (x_{1},x_{3})K_{F}(x_{3})\alpha _{4}\triangle x+0

=-{\frac  {1}{2}}K_{F}(x_{1})\alpha _{2}\triangle x+\sum _{{j=3}}^{4}\Gamma (x_{1},x_{{j-1}})K_{F}(x_{{j-1}})\alpha _{j}\triangle x

Similarly, for the third row, \pi ^{N}{\mathcal  {F}}(b^{N})(x_{2})=\gamma _{3}, and

\gamma _{3}={\frac  {1}{\triangle x}}\int _{{x_{2}}}^{{x_{3}}}\left[\int _{{x_{2}}}^{{x_{4}}}\Gamma (x_{2},y)K_{F}(y)b^{4}(y)dy\right]dx-{\frac  {1}{\triangle x}}\int _{{x_{2}}}^{{x_{3}}}{\frac  {1}{2}}K_{F}(x_{2})b^{4}(x_{2})\triangle xdx

=-{\frac  {1}{2}}K_{F}(x_{2})\alpha _{3}\triangle x+{\frac  {1}{\triangle x}}\int _{{x_{2}}}^{{x_{3}}}\left[\int _{{x_{2}}}^{{x_{3}}}\Gamma (x_{2},y)K_{F}(y)b^{4}(y)dy+\int _{{x_{3}}}^{{x_{4}}}\Gamma (x_{2},y)K_{F}(y)b^{4}(y)dy\right]dx

=-{\frac  {1}{2}}K_{F}(x_{2})\alpha _{3}\triangle x+{\frac  {1}{\triangle x}}\int _{{x_{2}}}^{{x_{3}}}\left[\Gamma (x_{2},x_{3})K_{F}(x_{3})b^{4}(x_{3})\triangle x+\Gamma (x_{2},x_{4})K_{F}(x_{4})b^{4}(x_{4})\triangle x\right]dx

=-{\frac  {1}{2}}K_{F}(x_{2})\alpha _{3}\triangle x+\Gamma (x_{2},x_{3})K_{F}(x_{3})\alpha _{4}\triangle x+0

Finally, for the fourth row, \pi ^{N}{\mathcal  {F}}(b^{N})(x_{3})=\gamma _{4}, and

\gamma _{4}={\frac  {1}{\triangle x}}\int _{{x_{3}}}^{{x_{4}}}\left[\int _{{x_{3}}}^{{x_{4}}}\Gamma (x_{3},y)K_{F}(y)b^{4}(y)dy\right]dx-{\frac  {1}{\triangle x}}\int _{{x_{3}}}^{{x_{4}}}{\frac  {1}{2}}K_{F}(x_{3})b^{4}(x_{3})\triangle xdx

=-{\frac  {1}{2}}K_{F}(x_{3})\alpha _{4}\triangle x+0

For the MATLAB program, we turned the compact column vector above, [\pi ^{N}{\mathcal  {F}}(b^{N})], into an n\times n matrix where the -{\frac  {1}{2}}K_{F}(x_{i})\triangle x terms sit nicely on the diagonal when we multiply by the density vector b=[\alpha _{1}\cdots \alpha _{N}]', and each term in the summations gets multiplied by its respective density term. The code to build \pi ^{N}{\mathcal  {F}}(b^{N}) follows.

% Build Projection of Fragmentation Matrix

% The logic here is to build a matrix 'Frag' such that Frag*b is % projection of Fragmentation as a function of the projected b

gammaF = 6.6e-5; % - parameter given in paper

K_F = @(xx) gammaF.*((xx - xmin).^(1/3));

Frag = zeros(N,N);

for m = 1:N

   Frag(m,m) = -.5*K_F(x(m));
   for k = m+1:N
       Frag(m,k) = probdist(x(m),x(k))*K_F(x(k));
   end

end

Frag(1,1) = 0; % - only difference in the pattern for the diagonal

% The whole matrix has a factor of delx

Frag = Frag.*delx;

Aggregation

Projected aggregation matrix, $\pi^N(\mathcal{A}(b^N))$

K_{A}(x,y) represents the aggregation kernel. For this analysis, we used the turbulent mixing kernel described by the authors as K_{A}(x,y)=\gamma _{A}\left({\frac  {\epsilon }{\nu }}\right)^{{1/2}}\left(x^{{1/3}}+y^{{1/3}}\right)^{3}, and we determine the following:

\pi ^{N}{\mathcal  {A}}(b^{N})]={\begin{bmatrix}-\alpha _{1}\sum _{{j=1}}^{{N-1}}K_{A}(x_{1},x_{j})\alpha _{j}\triangle x\\{\frac  {1}{2}}K_{A}(x_{1},x_{1})\alpha _{1}\alpha _{1}\triangle x-\alpha _{2}\sum _{{j=1}}^{{N-2}}K_{A}(x_{2},x_{j})\alpha _{j}\triangle x\\\vdots \\{\frac  {1}{2}}\sum _{{j=1}}^{{N-2}}K_{A}(x_{j},x_{{N-1-j}})\alpha _{j}\alpha _{{N-1-j}}\triangle x-\alpha _{{N-1}}K_{A}(x_{{N-1}},x_{1})\alpha _{1}\triangle x\\{\frac  {1}{2}}\sum _{{j=1}}^{{N-1}}K_{A}(x_{j},x_{{N-j}})\alpha _{j}\alpha _{{N-j}}\triangle x\end{bmatrix}}

For this matrix, we'll offer another perspective on how to understand the way in which we can derive it. In our discretization scheme, the first "bin" starts at x_{0}=2 fL and goes to x_{1}. We then take the approach that all aggregates in that bin have the discrete value infinitesimally close to x_{1}. Under that rationale, no two clumps can aggregate and produce one of size x_{0}, so in the first row, we can only reduce the density of size x_{0} aggregates, hence the strictly negative term. Furthermore, any size x_{0} group, which we represent in the matrix as size x_{1} (since it is infinitesimally close to size x_{1}) can aggregate with any other size group up to size x_{{N-2}} which we represent as size x_{{N-1}} (the size of the group at the large end of the x_{{N-2}} bin).

Another subtlety to highlight is that for each size x_{i} aggregate, we use the density for the x_{i} bin. In other words, in terms of size, we represent the x_{i} bin with size x_{{i+1}}, but we use the density of bin x_{i} which is \alpha _{i}. To illustrate, when an aggregate from the size x_{0} bin joins with an aggregate from the size x_{1} bin, we describe the interaction as \alpha _{1}K_{A}(x_{1},x_{2})\alpha _{2}, which creates an aggregate of not quite size x_{3} thereby increasing the density of aggregates of bin size x_{2}.

Finally, in our summation of the interactions between aggregates that increase the density of certain bin, we multiply by 1/2 to account for double counting. For instance, consider row three which accounts for the change in density of the size x_{2} bin. The terms that increase that density are represented by the sum,

{\frac  {1}{2}}\sum _{{j=1}}^{2}K_{A}(x_{j},x_{{3-j}})\alpha _{j}\alpha _{{3-j}}\triangle x={\frac  {1}{2}}\left(K_{A}(x_{1},x_{2})\alpha _{1}\alpha _{2}\triangle x+K_{A}(x_{2},x_{1})\alpha _{2}\alpha _{1}\triangle x\right)

The code to build \pi ^{N}{\mathcal  {A}}(b^{N}) follows.

% Build Projection of Aggregation Matrix

ee = 4.43e-4; % - Turbulent energy dissipation constant from paper

nu = 1.99e-6; % - Kinematic viscosity constant from paper

gammaA = 2.7e-15; % - parameter given in paper

% The sqrt of the ratio ee/nu is in 1/s, so I multiply by 60 to get in % 1/minutes

K_A = @(xx,yy) gammaA*sqrt(ee/nu)*60.*(xx.^(1/3) + yy.^(1/3)).^3;

Agg2 = zeros(N,1);

% This loop performs the second summation per my derivation

for m = 1:N-1;

   ka_b1 = K_A(x(m+1),x(2:N+1-m))'.*b(1:N-m);
   Agg2(m) = -b(m)*delx*sum(ka_b1);

end

% This loop performs the first summation per my derivation

for m = 2:N;

   for k = 2:m
       Agg2(m) = Agg2(m) + (.5*K_A(x(k),x(m-k+2))*b(k-1)*b(m-k+1)*delx);
   end

end

Results

Running our code from time, 0, to time, 240 minutes, we generate Figure 2.

Figure 2: Density Distribution with Aggregation, Fragmentation, and Growth.

We can also compare our results in Figure 3 with those of Figure 4 which were generated in the paper. Here we see similar shapes and orders of magnitude, but we must note that we are using a different initial density function.

Figure 3: Bin Distributions from this Wiki's Simulation.
Figure 4: Bin Distributions from the Paper.

In a case like this, where we do not necessarily have a final "true" answer that verifies whether we have the correct results, we can run some "sanity checks." For instance in Figure 5, we ran our code using only growth (no aggregation or fragmentation).

Figure 5: Density Distribution with Growth Only.

Here we set the initial condition to b_{0}(x_{0})=1 and b_{0}(x_{i})=0 for i=1,\dots ,N-1. Notice this simulation runs for 20000 minutes to verify that starting with a density of for bacteria of size 2 fL, we finish with with a density of one for bacteria of maximum size 1000 fL.

Conclusions

This Wiki sought to introduce the PDE model for the flocculation dynamics of bacteria in suspension that Bortz, et al. presented in their paper, "\textit{Klebsiella pneumoniae} Flocculation Dynamics", describe in detail a means for generating the necessary matrices for solving the discretized PDE, and most importantly generate a MATLAB code that actually solves the PDE. Based on our derivations in section 4, and our results along with our "sanity check" in section 5, we are confident we have met our goals.

Additionally, Bortz, et al., provide a "Fit Sensitivity" Analysis which helps indicate where improvements in the experiment can be made. While not addressed in this Wiki, we highly recommend reading this section for a solid example of the mathematical analysis behind design of experiments.

Appendix

A copy of the presentation for this wiki can be seen here, presentation .

Full Code:

close all; clear all; clc

%%%%%%%%%%%%%%%%%%%% Set Up

syms y

N = 50; % - number of steps from x_min to x_max %% If this changes, update

       % Flocdyn function as well

%Assumed min and max volumes with units (fL)

xmin = 2;

xmax = 1000;

delx = (xmax - xmin) / N; % - equispaced mesh size

x = [xmin:delx:xmax]; % - spatial discretization

%Inital density function with y representing continuous spatial variable

   % Using updated coefficients from Dr. Bortz (i.e. not the function in
   % the paper)

% - if this changes, the expression in intinitcond2 needs to change accordingly

b0 = 89.828*exp(-0.0035303*y) + 1946.6*1e6*exp(-1.3159*y);

% initalize density function after pi^N application from x_min to the last % step before x_max

b = zeros(1,N);

for s = 1:N;

   b(s) = intinitcond2(x(s),x(s+1))/delx;

end


%%%%%%%%%%%%% Solve ODE system % timespan

T = 240;

tspan = [0,T];

% call the solver

[t,B] = ode45(@Flocdyn4,tspan,b);

% plot the results

M = length(t);

set(axes,'ZScale','log')

set(gca,'XTick',[1:10:N-10+1,N]);

set(gca,'XTickLabel',[x(1):10*delx:2+((N-10)*delx),2+((N-1)*delx)]);

set(gca,'YTick',[1,M]);

set(gca,'YTickLabel',[0,T]);

view([65 30]);

title('Aggregation, Fragmentation, and Growth')

hold on

mesh(B)

xlabel('Volume (fL)');ylabel('Time (min)');zlabel('Density (number of particles/fL)');

% Find total number of aggregates to compare against Figure 5 in the paper % Initial time

binI(1) = sum(B(1,1:ceil(48/delx)))*delx - 3*B(1,1); % - vol's 5 to last vol less than 50

binI(2) = sum(B(1,ceil(48/delx)+1:ceil(98/delx)))*delx; % - vol's first above 50 to last below 100

binI(3) = sum(B(1,ceil(98/delx)+1:ceil(198/delx)))*delx; % - vol's 1st above 100 to last below 200

binI(4) = sum(B(1,ceil(198/delx)+1:ceil(298/delx)))*delx; % - vol's 1st above 200 to last below 300

binI(5) = sum(B(1,ceil(298/delx)+1:ceil(398/delx)))*delx; % - vol's 1st above 300 to last below 400

% Final time

binF(1) = sum(B(end,1:ceil(48/delx)))*delx - 3*B(end,1); % - vol's 5 to last vol less than 50

binF(2) = sum(B(end,ceil(48/delx)+1:ceil(98/delx)))*delx; % - vol's first above 50 to last below 100

binF(3) = sum(B(end,ceil(98/delx)+1:ceil(198/delx)))*delx; % - vol's 1st above 100 to last below 200

binF(4) = sum(B(end,ceil(198/delx)+1:ceil(298/delx)))*delx; % - vol's 1st above 200 to last below 300

binF(5) = sum(B(end,ceil(298/delx)+1:ceil(398/delx)))*delx; % - vol's 1st above 300 to last below 400

%Display the output

formatspec = 'Initial time:\t%0.5g\t%0.5g\t\t%0.5g\t\t%0.5g\t\t%0.5g\nFinal time:\t\t%0.5g\t%0.5g\t%0.5g\t%0.5g\t%0.5g';

str1 = sprintf('Bin Sizes:\t\t5-50\t\t50-100\t\t100-200\t\t200-300\t\t300-400');

str2 = sprintf(formatspec,binI,binF);

disp(str1)

disp(str2)

%plot these bins at initial and final times

bins = [50 100 200 300 400];

fig2 = figure

set(axes,'YScale','log')

set(gca,'XTick',bins);

set(gca,'XTickLabel','5-50|50-100|100-200|200-300|300-400');

hold on

plot(bins,binI,'--o',bins,binF,'-.x');

axis([0 450 10^3 10^8]); title('Distributions');

xlabel('Volume Ranges (in fL)'); ylabel('Total Number of Aggregates');

legend('Model Outputs Initial Time', 'Model Outputs Final Time');

Called Functions

  • ODE Solver

function [db_dt]= Flocdyn4(t,b)

N = 50; % - number of steps from x_min to x_max

%Assumed min and max volumes with units (fL)

xmin = 2;

xmax = 1000;

delx = (xmax - xmin) / N; % - equispaced mesh size

x = [xmin:delx:xmax]; % - spatial discretization


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Build Growth Matrix

gammaG = 6.8e-4; % - parameter given

G = @(xx) gammaG.*xx.*(1-(xx./xmax));

Growthvec = zeros(N,N);

Growthvec(1,1) = -G(x(2));

for r = 2:N

   Growthvec(r,r-1) = G(x(r));
   Growthvec(r,r) = - G(x(r+1));

end

% The whole matrix has a factor of 1/delx

Growthvec = Growthvec./delx;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Build Projection of Fragmentation Matrix

% The logic here is to build a matrix 'Frag' such that Frag*b is

% projection of Fragmentation as a function of the projected b

gammaF = 6.6e-5; % - parameter given in paper

K_F = @(xx) gammaF.*((xx - xmin).^(1/3));

Frag = zeros(N,N);

for m = 1:N

   Frag(m,m) = -.5*K_F(x(m));
   for k = m+1:N
       Frag(m,k) = probdist(x(m),x(k))*K_F(x(k));
   end

end

Frag(1,1) = 0; % - only difference in the pattern for the diagonal

% The whole matrix has a factor of delx

Frag = Frag.*delx;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Build Projection of Aggregation Matrix

ee = 4.43e-4; % - Turbulent energy dissipation constant from paper

nu = 1.99e-6; % - Kinematic viscosity constant from paper

gammaA = 2.7e-15; % - parameter given in paper

% The sqrt of the ratio ee/nu is in 1/s, so I multiply by 60 to get in % 1/minutes

K_A = @(xx,yy) gammaA*sqrt(ee/nu)*60.*(xx.^(1/3) + yy.^(1/3)).^3;

Agg2 = zeros(N,1);

% This loop performs the second summation per my derivation

for m = 1:N-1;

   ka_b1 = K_A(x(m+1),x(2:N+1-m))'.*b(1:N-m);
   Agg2(m) = -b(m)*delx*sum(ka_b1);

end

% This loop performs the first summation per my derivation

for m = 2:N;

   for k = 2:m
       Agg2(m) = Agg2(m) + (.5*K_A(x(k),x(m-k+2))*b(k-1)*b(m-k+1)*delx);
   end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Form the system of ODEs

db_dt = Growthvec*b + Agg2 + Frag*b;

  • Projecting Initial Conditions

function c = intinitcond2(x1,x2)

% This function will integrate the symbolic initial function from volume % x1 to volume x2.

c = 19466000000000/(13159*exp((13159*x1)/10000)) - ...

   19466000000000/(13159*exp((13159*x2)/10000)) + ...
   809098694654873829376/(31798115529012125*...
   exp((254384924232097*x1)/72057594037927936)) - ...
   809098694654873829376/(31798115529012125*...
   exp((254384924232097*x2)/72057594037927936)) ;

References

  1. D. M. Bortz, T. L. Jackson, K. A. Taylor, A. P. Thompson,and J. G. Younger, Klebsiella pneumoniae flocculation dynamics. Bulletin of Mathematical Biology, 70(3):745-768, April 2008.