February 19, 2017, Sunday
University of Colorado at Boulder Search A to Z Campus Map University of Colorado at Boulder CU 
Search Links

Dua Chaker

From MathBio

Jump to: navigation, search

Dua is doing research with Dr. Bortz

Bacteremia and Biofilms

Various types of bacteria exist and reside on the skin's surface and the face of internal organs that prove to be often harmless as they are sustained by the body's immune system. When the bacteria makes its way into the circulatory system, however, bacteremia symptoms usually arise. Bacteremia, the presence of bacteria in the blood stream, is more common in individuals with weaker immune systems often in cases of prolonged illness as well as those undergoing medical procedures. When the bacteria collectively adhere on the internal vessels in the blood stream, biofilms may develop. Through bacteria colonization, the biofilm, protected and held by a extracellular matrix, begins to grow and detach eventually dispersing and spreading to various areas in the body. These biolfilms thus prove harmful as the immune system's response deems ineffective as it fails to penetrate the biofilm surface.

The focus then for the research is to model biofilms, with actual bactermia locations in the blood, using mass spring systems. By representing the connections between bacteria with springs, various cases can be tested. By starting off with simple arbitrary scenarios in two dimensions and then analyzing the bacteria in three dimensions for more realistic cases, I was able to get an understanding of the behavior of the bacteria. Ultimately, the goal would be to then have a model for a biofilm that can be applied and tested for different initial conditions.

Modeling Bacteria in 3D

I analyzed various cases of bacteria in the blood stream by looking at different positions of the bacteria as well as the connection of the springs holding them together. By analyzing the damped mass spring system, I was able to look at several cases and come to different results. The first case, a single bacterium connected to the vessel wall by a spring at an angle to the right:

Figure 1. Case 1

By taking into account the sum of the forces on the bacteria in the x, y and z direction, and using the drag force equation, the spring force can be determined. In fact, I developed a function in Matlab that will determine the force acting on the spring given any arbitrary bacteria position and the position of the other end of the spring. This proves very helpful as this can be used for even more complicated cases that I will later discuss. (See Appendix for Matlab Code and mass spring system equations)

In addition to determining the spring force, I plotted the system and observed that the spring stretches in the positive x direction and also in the positive y direction. The velocity also increased in the positive direction. The second case I looked at, a single bacterium connected to the vessel wall by a spring at an angle to the left (keeping the origin at the initial clamp of the first case):

Figure 2. Case 2

Taking similar procedures as with the first case, like results followed as the spring also stretched out in the positive direction, but stretched out further. The velocity was also increased in the second case, the two cases compare as follows:

S1vsS2.jpg S1vsS22.jpg

After getting a general idea of the first 2 cases, I combined the two springs together and observed the behavior of the bacteria. The bacteria stretches out considerably more than the individual cases as there are more forces acting on the bacteria resulting in a larger stretch as well as higher velocity.

Case3vsSp12.jpg Case3vsSp122.jpg

Lastly, I looked at a case where instead of analyzing one bacterium, there were two bacteria connected by springs on both sides as well as a fifth spring connecting the two bacteria together:

Figure 3. Case 5

By looking at the sum of forces on each bacteria and creating equations for the mass spring systems, I was able to determine the exact forces acting on each spring and thus the x y and z components of each spring force. The following plots result:

Case4P.jpg Case4V.jpg

Both bacteria stretch in the positive direction, however even more in the case of the second bacteria. Keeping in mind that the laminar flow in the bloodstream is towards the right, both bacteria sway to the right, and since the second bacteria is further to the right than the first bacteria, it makes sense that the second bacteria stretches out longer as exhibited from the plots above.

Now that arbitrary cases of bacteria and spring positions have been analyzed and an idea of the behavior of bacteria in three dimensions has been developed, I take to applying the functions to real life cases. I begin by obtaining a set of data, representing three 2D sections of biofilm, which illustrates the exact bacterial positions in the blood stream. Given these sets of bacterial positions and utilizing the Delaunay Triangulation, the following figures results:

BacteriaPos1.jpg DelT1.jpg

BacteriaPos2.jpg DelT2.jpg

BacteriaPos3.jpg DelT3.jpg

The Delaunay Triangulation thus constructs connections between the points (bacteria positions) as noted from the graphs. By representing these connections as springs, the spring force for each spring can be determined depending on what two nodes make up the connection. In order to determine what points are connected and those that are not, an adjacency matrix can be utilized. An adjacency matrix is a useful tool that allows for a way to represent which vertices (nodes) of a plot (in this case triangulation) are adjacent to other vertices in the plot. The nxn matrix is essentially comprised of zeros and ones, with zeros along the diagonal. If there exists a connection between two points, a 1 is placed in the index of the matrix, but if there is no connection, it is denoted with a 0 in the matrix. For example, if there is a connection between points 1 and 2, there would be a 1 in the (1,2) and (2,1) index.

Keeping all this in mind, a function can be created that will compute the adjacency matrix given the triangulation (see Appendix for code) and once that is known, a function is made that will scan its corresponding adjacency matrix and determine what points are connected (see Appendix) and compute the spring force of the spring connection (see Appendix). If no such connection exists, a message will be displayed that no connection is made. This allows for a very useful tool, as given any 2D positions of biofilm in the blood, connections can be identified and the spring forces can be computed.

Thus, by being able to find the connections of a given set of bacterial positions and determining the spring forces of each connection, the biofilm can be simulated. Given all the forces acting on each bacteria, one can then simulate the biofilm and observe the behavior of the bacteria collectively as it exists in the blood stream. This would then be applied to any biofilm model and predictions as to what would happen to bacteremia can be studied. This would prove helpful as these generalizations of bacteria movement can pave the way for further research on how bacteremia can be treated and possible removed.


Matlab Code:

   function Fsp = sf(ba,s,Lrest)        %function that finds the spring force
    k = 2.2e-7; %spring constant (N/m)
    m = -k *(sqrt((s(1)- ba(1))^2+(s(2)- ba(2))^2+(s(3)- ba(3))^2)- Lrest);   %determines magnitude of spring force
   %b= vector position of the bacteria
   %s= vector position of where the spring is attached to the wall
     v = s - ba;   %direction of force
     w = v/sqrt(v(1)^2 + v(2)^2 + v(3)^2);      %normalizing v
     Fsp = m*w;    %vector force of spring

   function dxdt = bacteria3dSF(t,x)
   r = 2e-3; %mm
   r1 = 12.5; %radius of aorta = 12.5mm. radius of capillaries = .004 mm
   r3 = 1e-3; %mm = radius of bacteria
   m = 1e-12; %mass of bacteria (g)
   k = 2.2e-7; %spring constant (N/m)
   b = 1e-10; %damping constant (Ns/m)

   % r2 = r1 - yd;
   v_m = 4250; %mm/s maximum velocity. 4.25e3 mm/s http://www.rwc.uc.edu/koehler/biophys/3e.html
   d = 1e-3; %g/mm^3 density of fluid (water)
   % r1 = 0.125 mm. Large radius (of Aorta). http://www.answers.com/topic/blood-vessel
   % r2 = Distance from center of tube radius
   % r3 = 1e-1 mm. Radius of Bacteria. http://physics.bu.edu/~redner/211-sp06/class01/MKS.html
   Cd = 0.48; %drag coefficient for a sphere

   c = v_m/(r1^2); %Constant representing - 1/4n * deltaP/deltax

   %v = c*(r1^2 - r2^2); %Finding the velocity based on Poiseuille's equation
    %Fd = 0.5*d*(v^2)*Cd*(pi*r3^2);
   %s1 = [0 0 0]' ;
   %s2 = [4e-3 0 0]';
   %s3 = [4e-3 0 0]';
   %s4 = [8e-3 0 0]';
   ba = [1e-4 2e-3 0 3e-4 2e-3 0]';  %compiled position vector for 2 bacteria positions:- 1st bacteria position: [1e-4 2e-3 0]      
   and 2nd bacteria: [3e-4 2e-3 0]
   s = [0 0 0 4e-3 0 0 4e-3 0 0 8e-3 0 0]';   %compiled position vector for 4 spring ends positions:- 1st end position:[0 0 
   0],2nd end:[4e-3 0 0],3rd end: [4e-3 0 0],4th end:[ 8e-3 0 0]
   Lrest = 2e-3;    %rest length
   Fs1= sf([ba(1) ba(2) ba(3)]',[s(1) s(2) s(3)]',Lrest);  %spring force for spring #1 (acting on bacteria 1)
   Fs2= sf([ba(1) ba(2) ba(3)]',[s(4) s(5) s(6)]',Lrest);  %spring force for spring #2 (acting on bacteria 1)
   Fs3= sf([ba(4) ba(5) ba(6)]',[s(7) s(8) s(9)]',Lrest);  %spring force for spring #3 (acting on bacteria 2)
   Fs4= sf([ba(4) ba(5) ba(6)]',[s(10) s(11) s(12)]',Lrest);   %spring force for spring #4 (acting on bacteria 2)
   Fs5a= sf([ba(4) ba(5) ba(6)]',[ba(1) ba(2) ba(3)]', Lrest); %spring force for spring #5 (acting on bacteria 1)
   Fs5b= sf([ba(1) ba(2) ba(3)]',[ba(4) ba(5) ba(6)]', Lrest); %spring force for spring #5 (acting on bacteria 2)
   %bacteria 1 = x(1:6,1);
   %bacteria 2 = x(7:12,1);
     dxdt = zeros(12,1);                                                                     %system for bacteria 1
     dxdt(1) = x(2);
     dxdt(2) = (0.5*d*((c*(r1^2 - (r1 - x(3))^2))^2)*Cd*(pi*r3^2)+ Fs1(1) + Fs2(1)+ Fs5a(1)- b*x(2))/m;
     dxdt(3) = x(4);
     dxdt(4) = (Fs1(2) + Fs2(2)+ Fs5a(2) - b*x(4))/m;
     dxdt(5) = x(6);
     dxdt(6) = (Fs1(3) + Fs2(3)+ Fs5a(3) - b*x(6))/m;
     dxdt(7) = x(8);                                                                         %system for bacteria 2
     dxdt(8) = (0.5*d*((c*(r1^2 - (r1 - x(9))^2))^2)*Cd*(pi*r3^2)+ Fs3(1) + Fs4(1)+ Fs5b(1) - b*x(8))/m;
     dxdt(9) = x(10);
     dxdt(10) = (Fs3(2) + Fs4(2) + Fs5b(2) - b*x(10))/m;
     dxdt(11) = x(12);
     dxdt(12) = (Fs3(3) + Fs4(3) + Fs5b(3) - b*x(12))/m;

    [t,x]=ode15s(@bacteria3dSF, [0 .01], [1e-4 0 2e-3 0 0 0 3e-4 0 6e-3 0 0 0]);
    plot3(x(:, 1),x(:,3),x(:,5),'--k');   %plots bacteria position for Bacteria 1
    xlabel('x position  (mm)');
    ylabel('y position  (mm)');
    zlabel('z position  (mm))');
    title('aorta: 3d position Bacteria 1 and Bacteria 2, Case 4');
    hold on
    plot3(x(:, 7),x(:,9),x(:,11),'r');   %plots bacteria position for Bacteria 2
    legend('Bacteria 1','Bacteria 2');
    hold off

   plot3(x(:, 2),x(:,4),x(:,6),'--k');   %plots bacteria velocity for Bacteria 1
   xlabel('x velocity  (mm/s)');
   ylabel('y velocity  (mm/s)');
   zlabel('z velocity  (mm/s)');
   title('aorta: 3d velocity Bacteria 1 and Bacteria 2, Case 4');
   hold on
   plot3(x(:, 8),x(:,10),x(:,12),'r');   %plots bacteria velocity for Bacteria 2
   legend('Bacteria 1','Bacteria 2');
   hold off

Case 1 and 2 equations:

                                                          mx"=F_d+ F_(s1,i)  
                                                          my"= F_(s1,j) 
                                                          mz"= F_(s1,k)
                                                F_s1=-k*(√(x^2+y^2+z^2 )-L_0)

Case 3 equations:

                                                           mx"=F_d+ F_(s1,i)+ F_(s2,i)  
                                                           my"= F_(s1,j)+ F_(s2,j) 
                                                           mz"= F_(s1,k)+ F_(s2,k)
                                                F_s1=-k*(√(x^2+y^2+z^2 )-L_0)
                                                F_s2=-k*(√((x-x_d1 )^2+y^2+z^2 )-L_0)

Case 4 equations:

                                                           mx"=F_d+ F_(s1,i)+ F_(s2,i)+ F_(s5a,i)  
                                                           my"= F_(s1,j)+ F_(s2,j)+ F_(s5a,j)
                                                           mz"= F_(s1,k)+ F_(s2,k)+ F_(s5a,k)
                                                         mx_2 "=F_d+F_(s3,i)+ F_(s4,i)+ F_(s5b,i)
                                                         my_2 "= F_(s3,j)+ F_(s4,j)+ F_(s5b,j)
                                                         mz_2 "= F_(s3,k)+ F_(s4,k)+ F_(s5b,k)
                                                         F_s1=-k*(√(x^2+y^2+z^2 )-L_0)
                                                         F_s2=-k*(√((x-x_d1 )^2+y^2+z^2 )-L_0)
                                                         F_s3=-k*(√((x-x_d1 )^2+y^2+z^2 )-L_0)
                                                         F_s4=-k*(√((x-(x_d1+x_d2))^2+y^2+z^2 )-L_0)
                                                         F_s5=-k*(√((x-x_2))^2+(y-y_2))^2+(z-z_2))^2 )-L_0)

Delaunay Triangulation

       %Various plots for the data sets X,X2,X3
       %Three 2D sections of biofilm
       %Units are in microns, representing a 1 micron slice out of a 3D set of points
         plot(X(:,1),X(:,2),'bs','LineWidth',2,...    %plots the bacteria positions
       xlabel('x position');
       ylabel('y position');
       title('Set 1 - X');
       xlabel('x position');
       ylabel('y position');
       title('Set 1 - X');       
      TRI = delaunay(X);     %X= <130x2>    %delaunay triangulation for X
      dt = DelaunayTri(X);
      xlabel('x position');
      ylabel('y position');
      title('Set 1 - X');
       %hold on                           %labels the points and triangles
       %vxlabels = arrayfun(@(n) {sprintf('P%d', n)}, (1:130)');
       %Hpl = text(X(:,1),X(:,2),vxlabels, 'FontWeight', 'bold', 'HorizontalAlignment',...
       %'center', 'BackgroundColor', 'none');
       %ic = incenters(dt);
       %numtri = size(dt,1);
       %trilabels = arrayfun(@(X) {sprintf('T%d', X)}, (1:numtri)');
       %Htl = text(ic(:,1), ic(:,2), trilabels, 'FontWeight', 'bold', ...
       %'HorizontalAlignment', 'center', 'Color', 'blue');
       %hold off
        figure(4)            %representation plot in 3d
        title('Set 1 - X');
        disp('The triangulation datastructure for X is:')
        disp('The edges of the triangulation for X is:')
         plot(X2(:,1),X2(:,2),'bs','LineWidth',2,...   %plots the bacteria postions
         xlabel('x position');
         ylabel('y position');
         title('Set 2 - X2');
         TRI2 = delaunay(X2); %X2= <127x2>     %delaunay triangulation for X2
         dt2 = DelaunayTri(X2);
         xlabel('x position');
         ylabel('y position');
         title('Set 2 - X2');
          %hold on                              %labels the points and triangles
          %vxlabels2 = arrayfun(@(n) {sprintf('P%d', n)}, (1:127)');
          %Hp2 = text(X2(:,1),X2(:,2),vxlabels2, 'FontWeight', 'bold', 'HorizontalAlignment',...
          %'center', 'BackgroundColor', 'none');
          %ic2 = incenters(dt2);
          %numtri2 = size(dt2,1);
          %trilabels2 = arrayfun(@(X2) {sprintf('T%d', X2)}, (1:numtri)');
          %Ht2 = text(ic(:,1), ic(:,2), trilabels2, 'FontWeight', 'bold', ...
          %'HorizontalAlignment', 'center', 'Color', 'blue');
          %hold off
          %figure(7)                                %representation plot in 3d
           disp('The triangulation datastructure for X2 is:')
           disp('The edges of the triangulation for X2 is:')
           plot(X3(:,1),X3(:,2),'bs','LineWidth',2,...  %plots the bacteria postions
            xlabel('x position');
            ylabel('y position');
            title('Set 3 - X3');
            TRI3 = delaunay(X3);     %X3= <144x2>    %delaunay triangulation for X3
            dt3 = DelaunayTri(X3);
            xlabel('x position');
            ylabel('y position');
            title('Set 3 - X3');
           %hold on                                   %labels the points and triangles
           %vxlabels3 = arrayfun(@(n) {sprintf('P%d', n)}, (1:144)');
           %Hp3 = text(X(:,1),X(:,2),vxlabels, 'FontWeight', 'bold', 'HorizontalAlignment',...
           %'center', 'BackgroundColor', 'none');
           %ic3 = incenters(dt3);
           %numtri3 = size(dt3,1);
           %trilabels3 = arrayfun(@(X3) {sprintf('T%d', X3)}, (1:numtri)'); 
           %Ht3 = text(ic(:,1), ic(:,2), trilabels3, 'FontWeight', 'bold', ...
           %'HorizontalAlignment', 'center', 'Color', 'blue');
           %hold off
           figure(10)                                 %representation plot in 3d
           title('Set 3 - X3');
          disp('The triangulation datastructure for X3 is:')
          disp('The edges of the triangulation for X3 is:')

Adjacency Matrix

           %obtains the adjacency matrix from the delaunay triangulation 
            %adjacency matrix for bacterial positions 1 --> X = <130x2>
            nvert = max(max(dt.Triangulation));            
            nface = length(dt.Triangulation);
           A = zeros(nvert);
     for i=1:nface
           A(dt.Triangulation(i,1),dt.Triangulation(i,2)) = 1; A(dt.Triangulation(i,2),dt.Triangulation(i,3)) = 1;
           A(dt.Triangulation(i,3),dt.Triangulation(i,1)) = 1;
          % make sure that all edges are symmetric
           A(dt.Triangulation(i,2),dt.Triangulation(i,1)) = 1; A(dt.Triangulation(i,3),dt.Triangulation(i,2)) = 1;
           A(dt.Triangulation(i,1),dt.Triangulation(i,3)) = 1;
            %adjacency matrix for bacterial positions 2 --> X = <127x2>
            nvert2 = max(max(dt2.Triangulation));           
            nface2 = length(dt2.Triangulation);
            A2 = zeros(nvert2);
      for i=1:nface2
            A2(dt2.Triangulation(i,1),dt2.Triangulation(i,2)) = 1; A2(dt2.Triangulation(i,2),dt2.Triangulation(i,3)) = 1;
            A2(dt2.Triangulation(i,3),dt2.Triangulation(i,1)) = 1;
            % make sure that all edges are symmetric
           A2(dt2.Triangulation(i,2),dt2.Triangulation(i,1)) = 1; A2(dt2.Triangulation(i,3),dt2.Triangulation(i,2)) = 1;
           A2(dt2.Triangulation(i,1),dt2.Triangulation(i,3)) = 1;

           %adjacency matrix for bacterial positions 3 --> X = <144x2> 
           nvert3 = max(max(dt3.Triangulation));           
           nface3 = length(dt3.Triangulation);
           A3 = zeros(nvert3);
     for i=1:nface3
          A3(dt3.Triangulation(i,1),dt3.Triangulation(i,2)) = 1; A3(dt3.Triangulation(i,2),dt3.Triangulation(i,3)) = 1;
          A3(dt3.Triangulation(i,3),dt3.Triangulation(i,1)) = 1;
          % make sure that all edges are symmetric
          A3(dt3.Triangulation(i,2),dt3.Triangulation(i,1)) = 1; A3(dt3.Triangulation(i,3),dt3.Triangulation(i,2)) = 1;
          A3(dt3.Triangulation(i,1),dt3.Triangulation(i,3)) = 1;

Spring Forces

         %This will loop around the adjancency matrix of each set and determine
         %which points (nodes) are connected. In return, it will determine the spring
         %force of the spring between them. If certain points are not connected, a
         %message will be displayed.
  for j=1:130
   for i=1:130
  if A(i,j)|| A(j,i) ~= 0
     Conn1 = i;
     Conn2 = j;
     disp(['Connection between  ', num2str(Conn1), '  and  '  num2str(Conn2)])
     Lrest = 2e-3;
     Fspring = sfp(X(i,:),X(j,:),Lrest);
     disp(['spring force is ', num2str(Fspring)])
 elseif A(i,j)|| A(j,i) == 0
     NoC1 = i;
     NoC2 = j;
     disp([num2str(NoC1), '  and  '  num2str(NoC2), ' are not connected'])
 for j=1:127
   for i=1:127
  if A2(i,j)|| A2(j,i) ~= 0
     Conn11 = i;
     Conn22 = j;
     disp(['Connection between  ', num2str(Conn11), '  and  '  num2str(Conn22)])
     Lrest = 2e-3;
     Fspring2 = sfp(X2(i,:),X2(j,:),Lrest);
     disp(['spring force is ', num2str(Fspring2)])
 elseif A2(i,j)|| A2(j,i) == 0
     NoC11 = i;
     NoC22 = j;
     disp([num2str(NoC11), '  and  '  num2str(NoC22), ' are not connected'])
  for j=1:144
   for i=1:144
  if A3(i,j)|| A3(j,i) ~= 0
     Conn111 = i;
     Conn222 = j;
     disp(['Connection between  ', num2str(Conn111), '  and  '  num2str(Conn222)])
     Lrest = 2e-3;
     Fspring3 = sfp(X3(i,:),X3(j,:),Lrest);
     disp(['spring force is ', num2str(Fspring3)])
 elseif A3(i,j)|| A3(j,i) == 0
     NoC111 = i;
     NoC222 = j;
     disp([num2str(NoC111), '  and  '  num2str(NoC222), ' are not connected'])