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
- 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.
- Type of Model - This model attempts to obtain numerical rates for aggregation, fragmentation and growth of clusters of bacteria.
- 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
with and . Figure 1 shows a representation of the processes included in this model.
Figure 1: Aggregation, Fragmentation, and Growth Schematic.
Variables
- - volume of bacterial aggregate (in fL)
- , - assumed minimum and maximum volumes
- - time (in minutes unless otherwise stated)
- - density function representing the number of aggregates of size x per fL
- - assumed logistic growth rate of bacterial aggregates of size x
- - net rate at which bacterial flocs aggregate
- for
- for
- - assumed turbulent mixing kernel describing the rate at which flocs of volumes x and y join to make a floc of size x + y
- - net rate at which bacterial flocs fragment
- for
- for
- - describes rate at which flocs fragment
- - describes probability density of daughter flocs for the fragmentation of a parent floc of size y
Parameters
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, , 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 ,
where represents a node in the mesh from the minimum volume to the maximum volume . These functions form an orthogonal basis for
with projections
Finally, we need one more operator, , defined as
All of this provides the spatial discretization needed for MOL to provide the resulting system of ODEs
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, , 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.
The projection into of these initial conditions is derived as follows:
So we have,
Then to develop the left hand side of our system of differential equations recall that when . So for instance, , but for any . Now we have
In other words,
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 defined above for each calculation on the mesh. On a suggestion from Dr. Bortz, we ran a symbolic integration on using symbolic variables and to generate the expression for the code below. With a mesh of 50 increments from to , 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, , in our ODE system.
Approximating generators
of the infinitesimal generator
This was defined as
,
for . Recall, represents an assumed logistical growth with maximum size of 1000 fL, where . When the approximating growth generator is applied to the projected density function, , in our system of ODEs, we get the relation
To see why this is true, consider application of the approximating generator to , and note the substitution, . A key point to remember below is that from to and it is zero everywhere else. So when we integrate it across those nodes we get , whereas it integrates to zero across any other pair of nodes.
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
therefore
This shows that in row one of our vector for $\mathcal{G}^N b^N$ we have the negative contribution from of , and in row two, contributes . Let's consider the contributions from which should clearly illustrate the pattern for the rest of the vector for .
therefore
The pattern then emerges as follows,
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,
represents the fragmentation kernel. For this analysis, we used the kernel described by the authors as . Furthermore represents the authors' assumed uniform post-fragmetation distribution, centered around , with width,
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
To see how this matrix develops from the actual model, we will demonstrate an example where we set . Our partition then starts at fL and ends at fL. The fragmentation operator acting on is
and its projection is defined as
Now consider row one of , so
where,
For each of the inter integrals we'll use a right hand estimate. Also, note that and our basis elements for the prior node, i.e., otherwise it is zero. Then
Now consider row two where , and
Similarly, for the third row, , and
Finally, for the fourth row, , and
For the MATLAB program, we turned the compact column vector above, , into an matrix where the terms sit nicely on the diagonal when we multiply by the density vector , and each term in the summations gets multiplied by its respective density term. The code to build 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))$
represents the aggregation kernel. For this analysis, we used the turbulent mixing kernel described by the authors as , and we determine the following:
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 fL and goes to . We then take the approach that all aggregates in that bin have the discrete value infinitesimally close to . Under that rationale, no two clumps can aggregate and produce one of size , so in the first row, we can only reduce the density of size aggregates, hence the strictly negative term. Furthermore, any size group, which we represent in the matrix as size (since it is infinitesimally close to size ) can aggregate with any other size group up to size which we represent as size (the size of the group at the large end of the bin).
Another subtlety to highlight is that for each size aggregate, we use the density for the bin. In other words, in terms of size, we represent the bin with size , but we use the density of bin which is . To illustrate, when an aggregate from the size bin joins with an aggregate from the size bin, we describe the interaction as , which creates an aggregate of not quite size thereby increasing the density of aggregates of bin size .
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 bin. The terms that increase that density are represented by the sum,
The code to build 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 and for . 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
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
- ↑ 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.