September 20, 2014, Saturday
University of Colorado at Boulder Search A to Z Campus Map University of Colorado at Boulder CU 
Search Links

Kissler: Bacterial Floc Structure

From MathBio

Jump to: navigation, search

This page summarizes Stephen Kissler's work on characterizing bacterial floc structure.

Please note that this page is currently under construction.



This research focuses on determining the fractal dimension of bacterial flocs and characterizing floc density.

Euclidean Shapes, Natural Fractals, and Dimension

Figure 1 - Closely Packed Circles
Figure 1 - Closely Packed Circles
Figure 2 - Fractal Pattern
Figure 2 - Fractal Pattern

Many natural objects, ranging from mountain ranges, coastlines, chemical aggregates, and lightning strikes exhibit fractal behavior. Most of these 'natural fractals' do not display the same sort of exact self-similarity as shapes such as the Koch Snowflake, but instead have statistical self-similarity, in which some statistical measure of the object remains approximately the same on different scales. Bacterial flocs such as those treated in this summary exhibit this second type of self-similarity.

A second characteristic of fractals is a non-integer scaling power, sometimes referred to as a fractal dimension. If we examine the familiar case of a sphere, it is easily recognized that the volume of the sphere scales with the cube of its radius; mathematically, this may be written as V \sim r^3, where V is volume and r is radius. More generally, we may write the formula N=s^d, where N represents some quality of the object (previously Volume), s represents some characteristic length associated with the object (previously Radius), and d stands for the scaling power or dimension of N with respect to the characteristic length. We may easily solve for d, yielding d = \frac{log(N)}{log(s)}. Integer scaling powers (values of d) hold for Euclidean objects such as circles, spheres, cubes, etc. and for tightly-packed objects, such as the grid of circles in Figure 1. For this grid, note that N could correspond to the number of spheres in the grid, s could correspond to the length of one side of the grid in number of circles, and, if this were the case, d would be 2. Fractal objects, however, have scaling factors d that are non-integer values. For example, the fractal object in Figure 2 (borrowed from Environmental Transport Processes, 2009, P. 470) has a scaling factor or fractal dimension of 1.465. Note that this value is less than 2; it is characteristic of fractals that their fractal dimensions are smaller than what would be 'expected' for the dimension of a Euclidean or tightly-packed shape. We will examine this particular phenomenon with respect to true bacterial flocs.

Fractal Nature of Bacterial Flocs

Figure 3 - Fractal Dimension Plot for Floc 2
Figure 3 - Fractal Dimension Plot for Floc 2
Figure 4 - Fractal Dimension Plot for Floc 4
Figure 4 - Fractal Dimension Plot for Floc 4

Our experimental data consists of x,y,z-coordinates for 27 flocs of bacteria. The fractal dimension of this experimental data can be calculated from the slope of the best-fit line on a plot of the number of points inside spheres of various radii versus the radii of those spheres; see the middle portion of Figures 1 and 2 for examples of such plots. Let us first examine Figure 3; the left pane of this plot shows the plotted data for a bacterial floc. On the x-y plane lies a set of concentric circles, which are projections of the spheres used in the plot in the middle pane. Note that the axes in the middle pane display 'log(Number of Points)' vs 'log(Radius)'; the Radius term on the x-axis corresponds to radii of these spheres, and the Number of Points term corresponds to the number of bacteria within the sphere of that radius. All of the spheres are centered at the center of mass of the bacterial aggregate. In the middle pane, note the 'leveling out' of the data points as the radius becomes large. This happens when the spheres become so large that they no longer pick up new points; all of the points have already been included in a smaller sphere. Therefore, the best-fit line displayed on the middle pane does not take into account all of the points in the plot; rather, the user can choose which points to include in the linear regression. The right pane includes a residual plot for the linear fit, and the next paragraph will include more information on this data. These plots were generated using the fracdimCalc_2.m matlab file, included in the appendix below.

For most of the flocs, the residuals were sufficiently heteroscedastic to assume that a linear trend was a good fit for the data (for example, consider Figure 3). However, a few displayed some striking curved trends, such as Floc 4, included in Figure 4. Even though the line seems to fit the data well on the middle pane, a clear curve is apparent in the residuals plot, suggesting that the fractal dimension of the floc itself may be changing as we move outward from the center of the floc.

Floc Density

Figure 5 - Isosurface for Floc 4; Inverse Density Value = 5
Figure 5 - Isosurface for Floc 4; Inverse Density Value = 5

I used MATLAB's 'isosurface' tool to analyze the densities of the bacterial flocs. This tool does the equivalent of a contour plot, just with an extra dimension. First, I created a three-dimensional grid of points at regular intervals that formed a 'box' around the bacterial floc. From each of these points, I summed the inverse-distance from the point to each bacteria in the floc. That way, a short distance between the point and the bacteria would correspond to a high inverse-distance value and vice-versa. Using these inverse-distance values, I created for each floc a video showing the growth of the isosurface from the smallest inverse-distance value to the largest. An example of such a video may be seen at (include link); note that these videos show a progression from large to small isovalue, so the video first indicates the regions of highest density and progresses to show the regions of lowest density. Most of the flocs display high density near the center of mass, and the isosurface moves near-uniformly outward until it encloses the entire floc.




%Stephen Kissler

%The purpose of this code is to calculate the fractal dimension of a given
%set of points in three-space.  This is the second version of
%'fracdimCalcOriginal.m'.  A primary purpose of this version is to
%understand the steep slopes seen for some aggregates in

%The input 'data' is an nx3 column vector containing x,y,z coordinates of
%the points.

%Indeces used: E,

function[out1,out2] = fracdimCalc_2(points,numcircs)
%Count the number of points in the user input:
sizePoints = size(points);
numpoints = sizePoints(1);

%Plot the points:
hold on
xlabel('x position')
ylabel('y position')
zlabel('z position')

%Calculate center of mass, 'com':
xsum = sum(points(:,1));
ysum = sum(points(:,2));
zsum = sum(points(:,3));
xavg = xsum/numpoints;
yavg = ysum/numpoints;
zavg = zsum/numpoints;
com = [xavg,yavg,zavg]; %gives an ordered set for 'com', the center of mass
%note this assumes each bacteria has equivalent mass.

%Create pointsPolar array; first create a preliminary set of points,
%pointsPolarpre, that contains the input cartesian points, adjusted so that
%the center of mass is the origin:
pointsPolarpre = points;
pointsPolarpre(:,1) = pointsPolarpre(:,1) - com(1);
pointsPolarpre(:,2) = pointsPolarpre(:,2) - com(2);
pointsPolarpre(:,3) = pointsPolarpre(:,3) - com(3);

%Convert the points array into a set of polar data in [r,theta,phi]:
pointsPolar = zeros(length(points),3);
pointsPolar(:,1) = sqrt(pointsPolarpre(:,1).^2 + pointsPolarpre(:,2).^2 + ...
pointsPolar(:,2) = atan(pointsPolarpre(:,2)./pointsPolarpre(:,1)); %may have to be updated
%when it is possible to have x=0 and x,y<0.
pointsPolar(:,3) = acos(pointsPolarpre(:,3)./pointsPolar(:,1));

%Set to second subplot for the fractal dimension plot points:
hold on

%initialize array fracdimPts, an mx2 array where m is the number of radii
%taken into consideration; the first column contains the log of the current
%radius, and the second column contains the log of the number of points
%inside that radius:
fracdimPts = zeros(50,2);
%Find the max and min distances from the center of mass, to be used in later
maxRadius = max(pointsPolar(:,1));
minRadius = min(pointsPolar(:,1));

%fill in fracdimPts and plot; also plot radii on bacterial plot:
%Allow #numcircs radii (suggestion: 50) to be drawn, starting from the first radius at which a
%point is picked up and ending at the maximum radius.
for indexE = 1:numcircs
    %Create an nx1 array that indicates with a '1' which points are within the
    %given radius and a '0' those that are outside:
    indecesIncludedPoints = pointsPolar(:,1)<(minRadius+((indexE/numcircs)*(maxRadius-minRadius)));
    numIncludedPoints = sum(indecesIncludedPoints); %counts # of pts in given radius
    %the following plots a circle of the given radius on the bacterial plot:
    %Set bacterial plot as the current figure:
    %Set up a radius from 0 to 2*pi and a 'theta' array that contains
    %the given radius at each index, for the sake of plotting the circle:
    rngtheta = linspace(0,2*pi);
    rngradius = ones(1,length(rngtheta)).*(minRadius+((indexE/numcircs)*(maxRadius-minRadius)));
    %Convert to cartesian to shift the center of the circle to the center
    %of radius:
    xCartCirc = rngradius.*cos(rngtheta);
    yCartCirc = rngradius.*sin(rngtheta);
    %Shift center:
    xCartCirc = xCartCirc + com(1);
    yCartCirc = yCartCirc + com(2);
    %plot the points:
    hold on
    if indexE == 10
    elseif indexE == 20
    elseif indexE == 30
    elseif indexE == 40
    elseif indexE == 50
    %Return the current plot to the fractal dimension plot:
    %Create entry for fracdimPts:
    fracdimPts(indexE,:) = [log(minRadius+indexE/numcircs*(maxRadius-minRadius)),log(numIncludedPoints)];
    %Plot this point:
    %Number the point:
    hold on
%Allow the user to determine which points constitute a linear trend:
ptskeptLeft = input('Keep which points? start: ');
ptskeptRight = input('Keep which points? end: ');
fracdimPts = fracdimPts(ptskeptLeft:ptskeptRight,:); %discard 'leveled out points' from fracdimPts
fitcoeffs = polyfit(fracdimPts(:,1),fracdimPts(:,2),1); %Find linear fit of chosen points
slope = fitcoeffs(1); %picks out slope from the coefficients calculated in previous step
fitlinex = linspace(fracdimPts(1,1),fracdimPts(length(fracdimPts),1)); %for plotting
fitliney = fitcoeffs(1)*fitlinex + fitcoeffs(2); %for plotting; uses linear fit coefficients
%Plot the line whose slope corresponds to the fractal dimension:
ylabel('Log(#of points)')
title(['Fractional Dimension Plot: Slope = ',num2str(slope)])

%Set to third subplot to plot residuals for the linear fit:
%calculate residuals:
residuals = zeros(length(fracdimPts),1); %initialize variable
for indexG = 1:length(fracdimPts)
    %calculate each residual; store in array:
    residuals(indexG) = fracdimPts(indexG,2) - (fitcoeffs(1)*fracdimPts(indexG,1)+fitcoeffs(2));
%plot residuals
hold on
%plot x-axis for reference:
residx = linspace(fracdimPts(1,1),fracdimPts(length(fracdimPts),1));

%Adjust axes on bacterial plot:
axis([com(1)-1.3*maxRadius com(1)+1.3*maxRadius com(2)-1.3*maxRadius com(2)+1.3*maxRadius com(3)-1.3*maxRadius com(3)+1.3*maxRadius])

out1 = slope;
out2 = numpoints;


%Stephen Kissler

%This function finds the principal components of a data set, rotates the
%points so that these components lie along the x,y, and z axes, and creates
%a mesh ellipsoid around the points.  This is the evolution of PCA.m.

%Note that input data must be formatted in 3xn format, not the usual nx3 format.

%Indeces used: A,B,C

function[out1,out2] = PCAaligned(data)
%Set up preliminary variables:
X = data;
Xt = X';
numpoints = length(Xt);

%Center the data:
%Calculate center of mass, 'com':
xsum = sum(Xt(:,1));
ysum = sum(Xt(:,2));
zsum = sum(Xt(:,3));
xavg = xsum/numpoints;
yavg = ysum/numpoints;
zavg = zsum/numpoints;
com = [xavg,yavg,zavg]; %gives an ordered set for 'com', the center of mass
%note this assumes each bacteria has equivalent mass.

%Create 'centered' array; this contains the centered points:
centered = Xt;
centered(:,1) = centered(:,1) - com(1);
centered(:,2) = centered(:,2) - com(2);
centered(:,3) = centered(:,3) - com(3);

%Re-assign the data matrices:
Xt = centered;
X = centered';

%Find eigenvalues and eigenvectors of the centered points;
[V,D] = eig(X*Xt);

%Rotate the axes:
%Duplicate the eigenvectors in a variable we can work with:
evecs = V;
%Find theta to align the major axis with the x-axis, while still having z:
theta = acos([evecs(1,3) evecs(2,3) 0]*[1;0;0]/norm([evecs(1,3) evecs(2,3) 0]));
if evecs(2,3) < 0 %then allow for a clockwise rotation greater than pi radians:
    theta = 2*pi-theta;
%Rotate eigenvectors about the z-axis to account for this:
for indexC = 1:3
    evecs(:,indexC) = [cos(theta) sin(theta) 0;-sin(theta) cos(theta) 0; 0 0 1]*evecs(:,indexC);
%Find psi to align the major axis directly with the x-axis:
psi = acos([evecs(1,3) 0 evecs(3,3)]*[1;0;0]/norm([evecs(1,3) 0 evecs(3,3)]));
if evecs(3,3) > 0 %then allow for a clockwise rotation greater than pi radians:
    psi = 2*pi-psi;
%Rotate eigenvectors about the y-axis to account for this:
for indexD = 1:3
    evecs(:,indexD) = [cos(psi) 0 -sin(psi);0 1 0;sin(psi) 0 cos(psi)]*evecs(:,indexD);
%Find rho to align minor axis directly with z-axis:
rho = acos([0 evecs(2,1) evecs(3,1)]*[0;0;1]/norm([0 evecs(2,1) evecs(3,1)]));
if evecs(2,1) > 0 %then allow for a clockwise rotation greater than pi radians:
    rho = 2*pi-rho;
%Rotate eigenvectors about the x-axis to account for this:
for indexE = 1:3
    evecs(:,indexE) = [1 0 0;0 cos(rho) sin(rho);0 -sin(rho) cos(rho)]*evecs(:,indexE);

%Rotate the points:
%Perform rotation about z-axis for points:
for indexA = 1:length(X)
    X(:,indexA) = [cos(theta) sin(theta) 0;-sin(theta) cos(theta) 0; 0 0 1]*X(:,indexA);
%perform rotation about y-axis for points:
for indexB = 1:length(X)
    X(:,indexB) = [cos(psi) 0 -sin(psi);0 1 0;sin(psi) 0 cos(psi)]*X(:,indexB);
%perform rotation about x-axis for points:
for indexF = 1:length(X)
    X(:,indexF) = [1 0 0;0 cos(rho) sin(rho);0 -sin(rho) cos(rho)]*X(:,indexF);

%Redefine Xt with rotated points:
Xt = X';

%Calculate max displacement, the dist. of the furthest point from the origin:
radii = sqrt(Xt(:,1).^2 + Xt(:,2).^2 + Xt(:,3).^2);
maxradius = max(radii);

%Determine the factor by which axes must be scaled down; want the length of
%the longest axis to be equal to the max radius:
maxEigVal = max(max(D)); %find maximum eigenvalue
sf = sqrt(maxEigVal)/maxradius; %assign 'sf', or scaling factor.
%Now, note that maxEigVal/sf = maxradius.

%Plot the centered and rotated input data:
xlabel('X Position')
ylabel('Y Position')
zlabel('Z Position')
title('Bacterial Floc with PCA Axes')
hold on

%Plot the rotated principal axes:
plot3([0 evecs(1,1)*sqrt(D(1,1))/sf],[0 evecs(2,1)*sqrt(D(1,1))/sf],[0 evecs(3,1)*sqrt(D(1,1))/sf],'LineWidth',1)
hold on
plot3([0 evecs(1,2)*sqrt(D(2,2))/sf],[0 evecs(2,2)*sqrt(D(2,2))/sf],[0 evecs(3,2)*sqrt(D(2,2))/sf],'LineWidth',1.2)
plot3([0 evecs(1,3)*sqrt(D(3,3))/sf],[0 evecs(2,3)*sqrt(D(3,3))/sf],[0 evecs(3,3)*sqrt(D(3,3))/sf],'LineWidth',3)

plot3([0 -evecs(1,1)*sqrt(D(1,1))/sf],[0 -evecs(2,1)*sqrt(D(1,1))/sf],[0 -evecs(3,1)*sqrt(D(1,1))/sf],'LineWidth',1)
plot3([0 -evecs(1,2)*sqrt(D(2,2))/sf],[0 -evecs(2,2)*sqrt(D(2,2))/sf],[0 -evecs(3,2)*sqrt(D(2,2))/sf],'LineWidth',1.2)
plot3([0 -evecs(1,3)*sqrt(D(3,3))/sf],[0 -evecs(2,3)*sqrt(D(3,3))/sf],[0 -evecs(3,3)*sqrt(D(3,3))/sf],'LineWidth',3)
axis([-1.1*maxradius 1.1*maxradius -1.1*maxradius 1.1*maxradius -1.1*maxradius 1.1*maxradius])

%Give the nx3 array X', the coordinates of the rotated points, as the first
out1 = X';

%Give the 3x1 array of lengths of the principal axes, from major to minor,
%as the second output (note lenghts correspond to distance from origin to
%the tip of the axis):
out2 = abs([sqrt(D(3,3))/sf ; sqrt(D(2,2))/sf ; sqrt(D(1,1))/sf]);


%Stephen Kissler

%The purpose of this file is to use the isocube.m function to create a
%video of an isosurface with value expanding from high density to low
%density.  Output is a an array of movie frames (use MOVIE command to play)

function[out] = isoCubeVid(points)

%Calculate the densities at each point using the first part of isoCube.m:

%Run PCAaligned.m to align principal axes of the floc with the x,y,z
[rotatedPoints,axisLengths] = PCAaligned(points(:,1:3)');
%Note 'axisLengths' is a 3x1 array with the length of the major axis (from
%origin to tip of axis) in the (1,1) spot and the length of the minor axis
%(again from origin to tip of axis) in the (3,1) spot.

%Create the cube of points:
x = -axisLengths(1):5:axisLengths(1);
y = -axisLengths(2):5:axisLengths(2);
z = -axisLengths(3):5:axisLengths(3);
% % x = linspace(-axisLengths(1),axisLengths(1),100);
% % y = linspace(-axisLengths(2),axisLengths(2),100);
% % z = linspace(-axisLengths(3),axisLengths(3),100);

%Create meshgrid with x,y, and z to be used in isosurface plot:
[X,Y,Z] = meshgrid(x,y,z);

%Initialize a 3-dimensional array to contain the inverse distances at each
dists = zeros(length(y),length(x),length(z));

%Now, for each point, sum the inverse distances from the point to all
for indexA = 1:length(x)
    for indexB = 1:length(y)
        for indexC = 1:length(z)
            %find the sum of the inverse distances:
            %create an array to store distances to each point:
            InvDists = 1./sqrt((x(indexA)-rotatedPoints(:,1)).^2 + (y(indexB)-rotatedPoints(:,2)).^2 + (z(indexC)-rotatedPoints(:,3)).^2);
            sumInvDists = sum(InvDists);
            dists(indexB,indexA,indexC) = sumInvDists;

%Calculate the minimum and maximum densities:

minDensity = min(min(min(dists)));
maxDensity = max(max(max(dists)));

%Close all open figures (commenting out this command draws the isosurfae
%over the plot of the bacterial points and axes):
close all

%Create a video for 150 intermediate isovalues ranging from the min density
%to the max density, using the second part of isoCube.m:
for indexA = 1:500
    %create isosurface plot:
    p = patch(isosurface(X,Y,Z,dists,maxDensity-(maxDensity-minDensity)*indexA/500));
    isonormals(X,Y,Z,dists, p)
    set(p, 'FaceColor', 'blue', 'EdgeColor', 'none');
    axis([-1.1*axisLengths(1) 1.1*axisLengths(1) -1.1*axisLengths(1) 1.1*axisLengths(1) -1.1*axisLengths(1) 1.1*axisLengths(1)])
    xlabel('x position')
    ylabel('y position')
    zlabel('z position')
    camlight; lighting phong
    %Capture image:
    movFile(indexA) = getframe;
    %Delete plot to ready axes for the next isosurface:
out = movFile;