Difference between revisions of "Team:KU Leuven/Modeling/Hybrid"

Line 450: Line 450:
 
<div class="sidebarright">
 
<div class="sidebarright">
 
<div id="eq27text">
 
<div id="eq27text">
<p>
 
$X_1$: coordinate at next timestep [µm]<br/>
 
$x_1$: coordinate at current timestep [µm]<br/>
 
$\mu$: bacterial diffusion constant [µm^2/h]<br/>
 
$\Delta t$: timestep [h]<br/>
 
$N$: normal distribution [-]<br/>
 
</p>
 
</div>
 
</div>
 
 
<div class="sidebarright">
 
<div id="eq28text">
 
 
<p>
 
<p>
 
$X_1$: coordinate at next timestep [µm]<br/>
 
$X_1$: coordinate at next timestep [µm]<br/>
Line 1,797: Line 1,785:
 
which of course equals (22), corresponding to a scaled
 
which of course equals (22), corresponding to a scaled
 
version of the original intercellular distance.
 
version of the original intercellular distance.
</p>
 
<div id="eq28">
 
<p>
 
$$
 
\lambda = (\frac{x_1-x_2}{2 \cdot \sqrt{\mu \cdot \Delta t})^2
 
+ (\frac{y_1-y_2}{2 \cdot \sqrt{\mu \cdot \Delta t})^2 \\
 
\frac{R^2_0}{4 \cdot \mu \Delta t}
 
\;\;\; \text{(22)}
 
$$
 
</p>
 
</div>
 
<p>
 
 
We then
 
We then
 
define the standardized variable R’² as (eq. 21). To recap,
 
define the standardized variable R’² as (eq. 21). To recap,
Line 2,662: Line 2,638:
 
$("#eq26text").hide();
 
$("#eq26text").hide();
 
$("#eq27text").hide();
 
$("#eq27text").hide();
$("#eq28text").hide();
 
 
});
 
});
  
Line 2,698: Line 2,673:
 
$("#eq26text").hide();
 
$("#eq26text").hide();
 
$("#eq27text").hide();
 
$("#eq27text").hide();
$("#eq28text").hide();
 
 
});
 
});
  
Line 2,734: Line 2,708:
 
$("#eq26text").hide();
 
$("#eq26text").hide();
 
$("#eq27text").hide();
 
$("#eq27text").hide();
$("#eq28text").hide();
 
 
});
 
});
  
Line 2,770: Line 2,743:
 
$("#eq26text").hide();
 
$("#eq26text").hide();
 
$("#eq27text").hide();
 
$("#eq27text").hide();
$("#eq28text").hide();
 
 
});
 
});
  
Line 2,806: Line 2,778:
 
$("#eq26text").hide();
 
$("#eq26text").hide();
 
$("#eq27text").hide();
 
$("#eq27text").hide();
$("#eq28text").hide();
 
 
});
 
});
  
Line 2,842: Line 2,813:
 
$("#eq26text").hide();
 
$("#eq26text").hide();
 
$("#eq27text").hide();
 
$("#eq27text").hide();
$("#eq28text").hide();
 
 
});
 
});
  
Line 2,878: Line 2,848:
 
$("#eq26text").hide();
 
$("#eq26text").hide();
 
$("#eq27text").hide();
 
$("#eq27text").hide();
$("#eq28text").hide();
 
 
});
 
});
  
Line 2,914: Line 2,883:
 
$("#eq26text").hide();
 
$("#eq26text").hide();
 
$("#eq27text").hide();
 
$("#eq27text").hide();
$("#eq28text").hide();
 
 
});
 
});
  
Line 2,950: Line 2,918:
 
$("#eq26text").hide();
 
$("#eq26text").hide();
 
$("#eq27text").hide();
 
$("#eq27text").hide();
$("#eq28text").hide();
 
 
$(".sidebarleft").show();
 
$(".sidebarleft").show();
 
$("#eq9textleft").show();
 
$("#eq9textleft").show();
Line 2,993: Line 2,960:
 
$("#eq26text").hide();
 
$("#eq26text").hide();
 
$("#eq27text").hide();
 
$("#eq27text").hide();
$("#eq28text").hide();
 
 
});
 
});
  
Line 3,029: Line 2,995:
 
$("#eq26text").hide();
 
$("#eq26text").hide();
 
$("#eq27text").hide();
 
$("#eq27text").hide();
$("#eq28text").hide();
 
 
});
 
});
  
Line 3,065: Line 3,030:
 
$("#eq26text").hide();
 
$("#eq26text").hide();
 
$("#eq27text").hide();
 
$("#eq27text").hide();
$("#eq28text").hide();
 
 
});
 
});
  
Line 3,101: Line 3,065:
 
$("#eq26text").hide();
 
$("#eq26text").hide();
 
$("#eq27text").hide();
 
$("#eq27text").hide();
$("#eq28text").hide();
 
 
});
 
});
  
Line 3,137: Line 3,100:
 
$("#eq26text").hide();
 
$("#eq26text").hide();
 
$("#eq27text").hide();
 
$("#eq27text").hide();
$("#eq28text").hide();
 
 
});
 
});
  
Line 3,173: Line 3,135:
 
$("#eq26text").hide();
 
$("#eq26text").hide();
 
$("#eq27text").hide();
 
$("#eq27text").hide();
$("#eq28text").hide();
 
 
});
 
});
  
Line 3,209: Line 3,170:
 
$("#eq26text").hide();
 
$("#eq26text").hide();
 
$("#eq27text").hide();
 
$("#eq27text").hide();
$("#eq28text").hide();
 
 
$(".sidebarleft").show();
 
$(".sidebarleft").show();
 
$("#eq9textleft").hide();
 
$("#eq9textleft").hide();
Line 3,252: Line 3,212:
 
$("#eq26text").hide();
 
$("#eq26text").hide();
 
$("#eq27text").hide();
 
$("#eq27text").hide();
$("#eq28text").hide();
 
 
});
 
});
  
Line 3,288: Line 3,247:
 
$("#eq26text").hide();
 
$("#eq26text").hide();
 
$("#eq27text").hide();
 
$("#eq27text").hide();
$("#eq28text").hide();
 
 
});
 
});
  
Line 3,324: Line 3,282:
 
$("#eq26text").hide();
 
$("#eq26text").hide();
 
$("#eq27text").hide();
 
$("#eq27text").hide();
$("#eq28text").hide();
 
 
});
 
});
  
Line 3,360: Line 3,317:
 
$("#eq26text").hide();
 
$("#eq26text").hide();
 
$("#eq27text").hide();
 
$("#eq27text").hide();
$("#eq28text").hide();
 
 
});
 
});
  
Line 3,396: Line 3,352:
 
$("#eq26text").hide();
 
$("#eq26text").hide();
 
$("#eq27text").hide();
 
$("#eq27text").hide();
$("#eq28text").hide();
 
 
$(".sidebarleft").show();
 
$(".sidebarleft").show();
 
$("#eq9textleft").hide();
 
$("#eq9textleft").hide();
Line 3,439: Line 3,394:
 
$("#eq26text").hide();
 
$("#eq26text").hide();
 
$("#eq27text").hide();
 
$("#eq27text").hide();
$("#eq28text").hide();
 
 
$(".sidebarleft").show();
 
$(".sidebarleft").show();
 
$("#eq9textleft").hide();
 
$("#eq9textleft").hide();
Line 3,482: Line 3,436:
 
$("#eq26text").hide();
 
$("#eq26text").hide();
 
$("#eq27text").hide();
 
$("#eq27text").hide();
$("#eq28text").hide();
 
 
$(".sidebarleft").show();
 
$(".sidebarleft").show();
 
$("#eq9textleft").hide();
 
$("#eq9textleft").hide();
Line 3,525: Line 3,478:
 
$("#eq26text").hide();
 
$("#eq26text").hide();
 
$("#eq27text").hide();
 
$("#eq27text").hide();
$("#eq28text").hide();
 
 
});
 
});
  
Line 3,561: Line 3,513:
 
$("#eq26text").hide();
 
$("#eq26text").hide();
 
$("#eq27text").hide();
 
$("#eq27text").hide();
$("#eq28text").hide();
 
 
});
 
});
  
Line 3,597: Line 3,548:
 
$("#eq26text").show();
 
$("#eq26text").show();
 
$("#eq27text").hide();
 
$("#eq27text").hide();
$("#eq28text").hide();
 
 
});
 
});
  
Line 3,633: Line 3,583:
 
$("#eq26text").hide();
 
$("#eq26text").hide();
 
$("#eq27text").show();
 
$("#eq27text").show();
$("#eq28text").hide();
 
 
});
 
});
  
 
$( "#eq27" ).mouseleave(function() {
 
$( "#eq27" ).mouseleave(function() {
$(".sidebarright").hide();
 
});
 
 
$( "#eq28" ).mouseenter(function() {
 
$(".sidebarright").show();
 
$("#eq1text").hide();
 
$("#eq2text").hide();
 
$("#eq3text").hide();
 
$("#eq4text").hide();
 
$("#eq5text").hide();
 
$("#eq6text").hide();
 
$("#eq7text").hide();
 
$("#eq8text").hide();
 
$("#eq9text").hide();
 
$("#eq10text").hide();
 
$("#eq11text").hide();
 
$("#eq12text").hide();
 
$("#eq13text").hide();
 
$("#eq14text").hide();
 
$("#eq15text").hide();
 
$("#eq16text").hide();
 
$("#eq17text").hide();
 
$("#eq18text").hide();
 
$("#eq19text").hide();
 
$("#eq20text").hide();
 
$("#eq21text").hide();
 
$("#eq22text").hide();
 
$("#eq23text").hide();
 
$("#eq24text").hide();
 
$("#eq25text").hide();
 
$("#eq26text").hide();
 
$("#eq27text").hide();
 
$("#eq28text").show();
 
});
 
 
$( "#eq28" ).mouseleave(function() {
 
 
$(".sidebarright").hide();
 
$(".sidebarright").hide();
 
});
 
});
Line 3,679: Line 3,592:
 
</script>
 
</script>
 
</html>
 
</html>
{{KU_Leuven}}
 
{{KU_Leuven/css}}
 
{{KU_Leuven/Lightbox/css}}
 
<html>
 
 
    <!--load mathJax related stuff -->
 
    <script type="text/x-mathjax-config">
 
        MathJax.Hub.Config({tex2jax: {inlineMath: [['$','$'], ['\\(','\\)']]}});
 
        MathJax.Hub.Config({ SVG: { scale: 100 }});
 
        <!-- possible fonts TeX, STIX-Web, Asana-Math, Neo-Euler, Gyre-Pagella, Gyre-Termes and Latin-Modern. -->
 
        MathJax.Hub.Config({ SVG: { Font: "Asana-Math" }});
 
    </script>
 
    <script
 
        src="http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_SVG.js"
 
        type="text/javascript"></script>
 
 
    <script>
 
        $(document).onload(function () {
 
            $(".main-navm").hide();
 
        }
 
    </script>
 
 
<link rel="stylesheet" type="text/css"
 
href="https://2015.igem.org/Template:KU_Leuven/Lightbox/CSS?action=raw&ctype=text/css" />
 
<script type="text/javascript" src="https://2015.igem.org/Template:KU_Leuven/Javascript?&action=raw&ctype=text/javascript"></script>
 
 
    <style>
 
        #modeling {
 
            background-color: transparent;
 
            border-style: solid;
 
            border: 0 solid transparent;
 
            border-bottom: 5px solid #8b7a57;
 
        }
 
        #centernav:hover #modeling {
 
            /* this is active when your mouse moves is over the item */
 
            border: 0 solid transparent;
 
            border-bottom: 0 solid transparent;
 
        }
 
        .main-navm:hover #modeling {
 
            /* this is active when your mouse moves is over the item */
 
            border: 0 solid transparent;
 
            border-right: 0 solid transparent;
 
        }
 
        @media screen and (max-width: 1000px) {#modeling {
 
                border-bottom: 5px solid transparent;
 
                border-right: 5px solid #8b7a57;
 
            }
 
        }
 
        #content {
 
            background-color: transparent;
 
        }
 
 
.summaryimg{
 
opacity: 0.6;
 
}
 
    </style>
 
 
    <head>
 
        <link href="https://static.igem.org/mediawiki/2015/9/9c/Ku_Leuven_Favicon.gif"
 
            rel="icon"/>
 
        <link / </head
 
            href="https://static.igem.org/mediawiki/2015/9/9c/Ku_Leuven_Favicon.gif"
 
            rel="shortcut icon">
 
 
        <body>
 
<!-- sidebar texts -->
 
<div class="sidebarright">
 
<div id="eq1text">
 
<p>
 
$C$: concentration [nmol/cl]<br/>
 
$\vec{r}$: position vector [µm]<br/>
 
$t$: time [h]<br/>
 
$D$: diffusion constant [µm^2/h]<br/>
 
</p>
 
</div>
 
</div>
 
 
<div class="sidebarright">
 
<div id="eq2text">
 
<p>
 
$C$: concentration [nmol/cl]<br/>
 
$\vec{r}$: position vector [µm]<br/>
 
$t$: time [h]<br/>
 
$\alpha$: production constant [nmol/h]<br/>
 
$\rho_A$: density of type A bacteria [#/cl]<br/>
 
</p>
 
</div>
 
</div>
 
 
<div class="sidebarright">
 
<div id="eq3text">
 
<p>
 
$C$: concentration [nmol/cl]<br/>
 
$\vec{r}$: position vector [µm]<br/>
 
$t$: time [h]<br/>
 
$k$: degradation constant [1/h]<br/>
 
</p>
 
</div>
 
</div>
 
 
<div class="sidebarright">
 
<div id="eq4text">
 
<p>
 
$C$: concentration [nmol/cl]<br/>
 
$\vec{r}$: position vector [µm]<br/>
 
$t$: time [h]<br/>
 
$D$: diffusion constant [µm^2/h]<br/>
 
$\alpha$: production constant [nmol/h]<br/>
 
$\rho_A$: density of type A bacteria [#/cl]<br/>
 
$k$: degradation constant [1/h]<br/>
 
</p>
 
</div>
 
</div>
 
 
<div class="sidebarright">
 
<div id="eq5text">
 
<p>
 
$\vec{r}$: position vector [µm]<br/>
 
$t$: time [h]<br/>
 
$\vec{F}_{applied}$: applied force [mN]<br/>
 
$\gamma$: frictional coefficient [mN*h/µm]<br/>
 
</p>
 
</div>
 
</div>
 
 
<div class="sidebarright">
 
<div id="eq6text">
 
<p>
 
$\vec{r}$: position vector [µm]<br/>
 
$t$: time [h]<br/>
 
$\chi$: chemotactic sensitivity [µm^2*cl/(h*nmol)]<br/>
 
$S$: chemoattractant concentration [nmol/cl]<br/>
 
$\mu$: bacterial diffusion constant [µm^2/h]<br/>
 
$\vec{W}$: Wiener process [h^(1/2)]<br/>
 
$\kappa$: chemotactic sensitivity constant [-]<br/>
 
</p>
 
</div>
 
</div>
 
 
<div class="sidebarright">
 
<div id="eq7text">
 
<p>
 
$n$: bacteria density [#/cl]<br/>
 
$\vec{r}$: position vector [µm]<br/>
 
$t$: time [h]<br/>
 
$\mu$: bacterial diffusion constant [µm^2/h]<br/>
 
$\chi$: chemotactic sensitivity [µm^2*cl/(h*nmol)]<br/>
 
$S$: chemoattractant concentration [nmol/cl]<br/>
 
$\kappa$: chemotactic sensitivity coefficient [-]<br/>
 
</p>
 
</div>
 
</div>
 
 
<div class="sidebarright">
 
<div id="eq8text">
 
<p>
 
$\chi$: chemotactic sensitivity [µm^2*cl/(h*nmol)]<br/>
 
$S$: chemoattractant concentration [nmol/cl]<br/>
 
$\mu$: bacterial diffusion constant [µm^2/h]<br/>
 
$\kappa$: chemotactic sensitivity constant [-]<br/>
 
$S$: chemoattractant concentration [nmol/cl]<br/>
 
$\vec{r}$: position vector [µm]<br/>
 
$t$: time [h]<br/>
 
</p>
 
</div>
 
</div>
 
 
<div class="sidebarright">
 
<div id="eq9text">
 
<p>
 
$E_p$: interaction potential [nJ]<br/>
 
$r_{ij}$: intercellular distance [µm]<br/>
 
$k$: “spring constant” [mN/µm]<br/>
 
$R_{cutoff}$: cutoff intercellular distance [µm]<br/>
 
$R_0$: equilibrium intercellular distance [µm]<br/>
 
$C$: constant to ensure continuity [nJ]<br/>
 
</p>
 
</div>
 
</div>
 
 
<div class="sidebarleft">
 
<div id="eq9textleft">
 
<p>
 
$C_1=-\frac{1}{2}\cdot k_3 \cdot(R_{cutoff}-R_0)^2$ <br/>
 
$C_2=-\frac{1}{2}\cdot k_1 \cdot(\frac{R_0}{2}-\frac{k_1+k_2}{k_1}\cdot
 
\frac{R_0}{2})^2 \\
 
+\frac{1}{2}\cdot k_2 \cdot(\frac{R_0}{2}-R_0)^2
 
+C_1$<br/>
 
</p>
 
</div>
 
</div>
 
 
<div class="sidebarright">
 
<div id="eq10text">
 
<p>
 
$\vec{F}_i$: cell-cell force acting on cell i [mN]<br/>
 
$E_p$: interaction potential [nJ]<br/>
 
$r_{ij}$: intercellular distance [µm]<br/>
 
$\vec{e}_{ij}$: unit vector from cell i to cell j [-]<br/>
 
$k$: “spring constant” [mN/µm]<br/>
 
$R_{cutoff}$: cutoff intercellular distance [µm]<br/>
 
$R_0$: equilibrium intercellular distance [µm]<br/>
 
</p>
 
</div>
 
</div>
 
 
<div class="sidebarright">
 
<div id="eq11text">
 
<p>
 
$\vec{r}$: position vector [µm]<br/>
 
$t$: time [h]<br/>
 
$\mu$: bacterial diffusion constant [µm^2/h]<br/>
 
$H$: concentration of AHL [nmol/cl]<br/>
 
$\vec{W}$: Wiener process [h^(1/2)]<br/>
 
$\gamma$: frictional coefficient [mN*h/µm]<br/>
 
$E_p$: interaction potential [nJ]<br/>
 
$r_{ij}$: intercellular distance [µm]<br/>
 
$\vec{e}_{ij}$: unit vector from cell i to cell j<br/>
 
$\chi$: chemotactic sensitivity [µm^2*cl/(h*nmol)]<br/>
 
$L$: concentration of leucine [nmol/cl]<br/>
 
</p>
 
</div>
 
</div>
 
 
<div class="sidebarright">
 
<div id="eq12text">
 
<p>
 
$\chi$: chemotactic sensitivity [µm^2*cl/(h*nmol)]<br/>
 
$L$: concentration of leucine [nmol/cl]<br/>
 
$H$: concentration of AHL [nmol/cl]<br/>
 
$\mu$: bacterial diffusion constant [µm^2/h]<br/>
 
$\kappa$: chemotactic sensitivity constant [-]<br/>
 
$\vec{r}$: position vector [µm]<br/>
 
$t$: time [h]<br/>
 
</p>
 
</div>
 
</div>
 
 
<div class="sidebarright">
 
<div id="eq13text">
 
<p>
 
$\mu$: bacterial diffusion constant [µm^2/h]<br/>
 
$H$: concentration of AHL [nmol/cl]<br/>
 
$\vec{r}$: position vector [µm]<br/>
 
$t$: time [h]<br/>
 
</p>
 
</div>
 
</div>
 
 
<div class="sidebarright">
 
<div id="eq14text">
 
<p>
 
$K$: probability density [-]<br/>
 
$x$: stochastic variable [-]<br/>
 
</p>
 
</div>
 
</div>
 
 
<div class="sidebarright">
 
<div id="eq15text">
 
<p>
 
$K_h$: kernel density [#/µm]<br/>
 
$x$: position [µm]<br/>
 
$h$: bandwidth [µm]<br/>
 
</p>
 
</div>
 
</div>
 
 
<div class="sidebarright">
 
<div id="eq16text">
 
<p>
 
$\rho$: agent density [#/cl]<br/>
 
$\beta$: conversion factor [µm/cl]<br/>
 
$N$: number of agents [#]<br/>
 
$K_h$: kernel density [#/µm]<br/>
 
$x$: position [µm]<br/>
 
$x_i$: position of ith agent [µm]<br/>
 
</p>
 
</div>
 
</div>
 
 
<div class="sidebarleft">
 
<div id="eq16textleft">
 
<p>
 
Conversion factor $\beta$ is needed to express the 1-D density [#/µm]
 
in units used in the PDE model [#/cl].
 
For the 2-D KDE $\beta$ has units [µm^2/cl]. The value of $\beta$
 
depends on assumptions made about the dimensions of the physical domain
 
which are not modeled. For example, when modeling a 2-D virtual domain
 
this constant will be high if bacteria are densely distributed
 
along the height of the experimental petri dish.
 
<br/>
 
</p>
 
</div>
 
</div>
 
 
<div class="sidebarright">
 
<div id="eq17text">
 
<p>
 
$f$: function, e.g. concentration [nmol/cl]<br/>
 
$\hat{f}$: interpolated function, e.g. [nmol/cl]<br/>
 
$x$: independent variable, e.g. position [µm]<br/>
 
$x_i$: independent variable for ith measurement, e.g. [µm]<br/>
 
</p>
 
</div>
 
</div>
 
 
<div class="sidebarright">
 
<div id="eq18text">
 
<p>
 
$f$: function, e.g. concentration [nmol/cl]<br/>
 
$\hat{f}$: interpolated function, e.g. [nmol/cl]<br/>
 
$x$: independent variable, e.g. position [µm]<br/>
 
$x_i$: independent variable for measurement (i,j), e.g. [µm]<br/>
 
$y$: independent variable, e.g. position [µm]<br/>
 
$y_j$: independent variable for measurement (i,j), e.g. [µm]<br/>
 
$Q_{i,j}$: symbol for grid point corresponding to measurement (i,j)<br/>
 
</p>
 
</div>
 
</div>
 
 
<div class="sidebarright">
 
<div id="eq19text">
 
<p>
 
$f$: function, e.g. concentration [nmol/cl]<br/>
 
$\hat{f}$: interpolated function, e.g. [nmol/cl]<br/>
 
$x$: independent variable, e.g. position [µm]<br/>
 
$x_i$: independent variable for measurement (i,j), e.g. [µm]<br/>
 
$y$: independent variable, e.g. position [µm]<br/>
 
$y_j$: independent variable for measurement (i,j), e.g. [µm]<br/>
 
$Q_{i,j}$: grid point corresponding to measurement (i,j) [-]<br/>
 
</p>
 
</div>
 
</div>
 
 
<div class="sidebarright">
 
<div id="eq20text">
 
<p>
 
$Re$: Reynolds number [-]<br/>
 
$\rho$: mass density [kg/m^3]<br/>
 
$v$: speed [m/s]<br/>
 
$L$: characteristic length [m]<br/>
 
$\eta$: dynamic viscosity [Pa*s]<br/>
 
</p>
 
</div>
 
</div>
 
 
<div class="sidebarright">
 
<div id="eq21text">
 
<p>
 
$F$: force [N]<br/>
 
$m$: mass [kg]<br/>
 
$x$: position [m]<br/>
 
$t$: time [s]<br/>
 
$\eta$: dynamic viscosity [Pa*s]<br/>
 
$a$: radius of particle [m]<br/>
 
$v$: speed [m/s]<br/>
 
</p>
 
</div>
 
</div>
 
 
<div class="sidebarleft">
 
<div id="eq21textleft">
 
<p>
 
Stokes' law states that the force of viscosity for a
 
spherical object moving through a viscous fluid is given
 
by $F=6\pi\eta av$ with $a$ the radius of the sphere. Importantly,
 
this force is proportional to the velocity $v$.
 
<br/>
 
</p>
 
</div>
 
</div>
 
 
<div class="sidebarright">
 
<div id="eq22text">
 
<p>
 
$v$: speed [m/s]<br/>
 
$v_0$: initial speed [m/s]<br/>
 
$\eta$: dynamic viscosity [Pa*s]<br/>
 
$a$: radius of cell [m]<br/>
 
$m$: mass [kg]<br/>
 
$t$: time [s]<br/>
 
$t_0$: start time [s]<br/>
 
</p>
 
</div>
 
</div>
 
 
<div class="sidebarleft">
 
<div id="eq22textleft">
 
<p>
 
To calculate the mass of a bacterium, we assume that the
 
bacterium is spherical with radius $a$ and the density
 
is twice that of water. Therefore, $m=4/3 \cdot \pi \cdot a^3 \cdot (2 \rho)$.
 
</p>
 
</div>
 
</div>
 
 
<div class="sidebarright">
 
<div id="eq23text">
 
<p>
 
$\mu$: bacterial diffusion constant [µm^2/h]<br/>
 
$\vec{W}$: Wiener process [h^(1/2)]<br/>
 
$\Delta t$: timestep [h]<br/>
 
$\vec{N}$: normal distribution [-]<br/>
 
</p>
 
</div>
 
</div>
 
 
<div class="sidebarleft">
 
<div id="eq23textleft">
 
<p>
 
$\vec{N}$ indicates a vector of which the components
 
are independent random variables that are normally
 
distributed.
 
</p>
 
</div>
 
</div>
 
 
<div class="sidebarright">
 
<div id="eq24text">
 
<p>
 
$\mu$: bacterial diffusion constant [µm^2/h]<br/>
 
$\Delta t$: timestep [h]<br/>
 
$\vec{N}$: normal distribution [-]<br/>
 
</p>
 
</div>
 
</div>
 
 
<div class="sidebarright">
 
<div id="eq25text">
 
<p>
 
$X_1$: coordinate at next timestep [µm]<br/>
 
$x_1$: coordinate at current timestep [µm]<br/>
 
$\mu$: bacterial diffusion constant [µm^2/h]<br/>
 
$\Delta t$: timestep [h]<br/>
 
$N$: normal distribution [-]<br/>
 
</p>
 
</div>
 
</div>
 
 
<div class="sidebarright">
 
<div id="eq26text">
 
<p>
 
$R$: intercellular distance at next timestep [µm]<br/>
 
$X_1$: coordinate at next timestep [µm]<br/>
 
</p>
 
</div>
 
</div>
 
 
<div class="sidebarright">
 
<div id="eq27text">
 
<p>
 
$X_1$: coordinate at next timestep [µm]<br/>
 
$x_1$: coordinate at current timestep [µm]<br/>
 
$\mu$: bacterial diffusion constant [µm^2/h]<br/>
 
$\Delta t$: timestep [h]<br/>
 
$N$: normal distribution [-]<br/>
 
</p>
 
</div>
 
</div>
 
 
<!--Begin Content -->
 
            <div class="summaryheader">
 
                <div class="summaryimg">
 
                    <img src="https://static.igem.org/mediawiki/2015/5/5c/KU_Leuven_Banner_Groen2.jpg"
 
width="100%"></img>
 
                    <div class="head">
 
                        <h2>
 
                            The Hybrid Model
 
                        </h2>
 
</div><!-- head -->
 
</div><!-- summaryimg -->
 
</div><!-- summary header -->
 
 
            <div class="summarytext1">
 
                <div class="part">
 
<h2>Introduction</h2>
 
                    <p> The hybrid model represents an intermediate level of detail in between the
 
                        colony level model and the internal model. Bacteria are treated as individual agents
 
                        that behave according to Keller-Segel type stochastic differential
 
                        equations, while chemical species are modeled using partial differential
 
equations. These different models are implemented and coupled within a single
 
hybrid modeling framework.
 
</p>
 
</div>
 
 
                <div class="part">
 
                    <h2>Partial Differential Equations</h2>
 
                    <p>
 
                        Spatial reaction-diffusion models that rely on
 
                        Partial Differential Equations (PDEs) are based on the assumption that the
 
                        collective behavior of individual entities, such as molecules or bacteria, can
 
                        be abstracted to the behavior of a continuous field that represents the density
 
                        of those entities. The brownian motion of molecules, for instance, is the result
 
                        of inherently stochastic processes that take place at the individual molecule
 
                        level, but is modeled at the density level by Fick’s laws of diffusion. These
 
                        PDE-based models provide a robust method to predict the evolution of large-scale
 
                        systems, but are only valid when the spatiotemporal scale is sufficiently large
 
                        to neglect small-scale stochastic fluctuations and physical granularity.
 
                        Moreover, such a continuous field approximation can only be made if the behavior
 
                        of the individual entities is well described.
 
                    </p>
 
</div>
 
 
                <div class="part">
 
                    <h2>Agent-Based Models</h2>
 
                    <p>
 
                        Agent-based models on the other hand explicitly treat the entities as individual
 
                        “agents” that behave according to a set of “agent rules”. An agent is an object
 
                        that acts independently from other agents and is influenced only by its local
 
                        environment. The goal in agent-based models is to study the emergent
 
                        systems-level properties of a collection of individual agents that follow
 
                        relatively simple rules. In smoothed particle hydrodynamics for example, fluids
 
                        are simulated by calculating the trajectory of each individual fluid particle at
 
                        every timestep. Fluid properties such as the momentum at a certain point can
 
                        then be sampled by taking a weighted sum of the momenta of the surrounding fluid
 
                        particles. A large advantage of agent-based models is that the agent rules are
 
                        arbitrarily complex and thus they allow us to model systems that do not
 
                        correspond to any existing or easily derivable PDE model. However, because every
 
                        agent is stored in memory and needs to be processed individually, simulating
 
                        agent-based models can be computationally intensive.
 
                    </p>
 
</div>
 
 
                <div class="part">
 
                    <h2>Hybrid Modeling Framework</h2>
 
                    <p>
 
                      In our system, there are both bacteria and chemical species that spread out
 
                        and interact on a petri dish to form patterns. On the one hand, the bacteria are
 
                        rather complex entities that move along chemical gradients and interact with one
 
                        another. Therefore they are ideally modeled using an agent-based model. On the
 
                        other hand, the diffusion and dynamics of the chemicals leucine and AHL are
 
                        easily described by well-established PDEs. To make use of the advantages of each
 
                        modeling approach, we decided to combine these two different types of modeling
 
                        in a hybrid modeling framework. In this framework we modeled the bacteria as
 
                        agents, while the chemical species were modeled using PDEs. There were two
 
                        challenges to our hybrid approach, namely coupling the models and matching them.
 
                        Coupling refers to the transfer of information between the models and matching
 
                        refers to dealing with different spatial and temporal scales to achieve
 
                        accurate, yet computationally tractable simulations.
 
 
                        <br/>
 
                        <br/>
 
 
                        In the following paragraphs we first introduce our hybrid model and its
 
                        coupling. Once the basic framework is established, the agent-based module and
 
                        PDE module are discussed in more depth and the issue of matching is highlighted.
 
                        We also expand on important aspects of the model and its implementation such as
 
                        boundary conditions and choice of timesteps. Then the results for the 1-D model
 
                        and 2-D model simulations are shown and summarized. Finally, the incorporation
 
                        of the internal model into the hybrid model is discussed and a proof of concept
 
                        is demonstrated.
 
                    </p>
 
</div><!-- part -->
 
 
<div class="whiterow"></div>
 
          </div><!-- end of summarytext1 -->
 
 
          <div class="summaryheader">
 
                <div class="summaryimg">
 
                    <img src="https://static.igem.org/mediawiki/2015/5/5c/KU_Leuven_Banner_Groen2.jpg"
 
                        width="100%">
 
                    <div class="head">
 
                        <h2>
 
                          Model Description
 
                        </h2>
 
</div><!-- head -->
 
</div><!-- summaryimg -->
 
</div><!-- summaryheader -->
 
            <div class="summarytext1">
 
 
                <div class="part">
 
<h2>System</h2>
 
                    <p>The main protagonists in our pattern-forming system are cell types A and B,
 
                        AHL and leucine. Cells A produce AHL as well as leucine. They are unaffected by
 
                        leucine, while cells B are repelled by leucine. AHL modulates the motility of
 
                        both cell types A and B, but in opposite ways. High concentrations of AHL will
 
                        render cell type A unable to swim but will activate cell type B’s motility.
 
                        Conversely, low concentrations of AHL causes swimming of cell type A and
 
                        incessant tumbling (thus immobility) of cell type B. Lastly, cells A express the
 
                        adhesin membrane protein, which causes them to stick to each other. Simply said,
 
                        our system should produce spots of immobile, sticky groups of A type cells,
 
                        surrounded by rings of B type cells. Any cell that finds itself outside of the
 
                        region that it should be in, is able to swim to their correct place and becomes
 
                        immobile there. More details can be found in the
 
                        <a href="https://2015.igem.org/Team:KU_Leuven/Research">research section</a>.</p>
 
</div><!-- part -->
 
 
                <div class="part">
 
<h2>Partial Differential Equations</h2>
 
                    <p>
 
                        As discussed in the previous paragraph, our hybrid model incorporates chemical
 
                        species using PDEs. In our system these are AHL and leucine. The diffusion of
 
                        AHL and leucine can be described by the heat equation (1).
 
</p>
 
<div id="eq1">
 
<p>
 
$$\frac{\partial C(\vec{r},t)}{\partial t}=D \cdot \nabla^2 C(\vec{r},t) \;\;\;
 
\text{(1)}$$
 
</p>
 
</div>
 
<p>
 
                        By using (1)
 
                        we assume that the diffusion speed is isotropic, i.e. the same in all spatial
 
                        directions. This also explains why it is called the heat equation, since heat
 
                        diffuses equally fast in all directions. A detailed explanation of the heat
 
                        equation can be found in box 1. The second factor that needs be taken into
 
                        account is the production of AHL and leucine by type A bacteria. In principle,
 
                        AHL and leucine production is dependent on the dynamically-evolving internal
 
                        states of all cells of type A. However, for our hybrid model we ignored the
 
                        inner life of all bacteria and instead assumed that AHL and leucine production
 
                        is directly proportional to the density of A type cells (2).
 
</p>
 
<div id="eq2">
 
<p>
 
                        $$  \frac{\partial C(\vec{r},t)}{\partial t}=\alpha \cdot \rho_A(\vec{r},t)
 
                        \;\;\; \text{(2)}$$
 
</p>
 
</div>
 
<p>
 
                        In the last
 
                        paragraph we will reconsider this assumption and assign each cell an internal
 
model. Finally, AHL and leucine are organic molecules and therefore they will
 
degrade over time.
 
 
                        We assume first-order kinetics meaning
 
                        that the rate at which AHL and similarly leucine disappear is proportional to
 
their respective concentrations (3a and 3b) assuming neutral pH
 
<sup><a href="#Schaefer2000">[6]</a></sup>.
 
</p>
 
<div id="eq3">
 
<p>
 
                        $$ \frac{\partial C_{AHL}(\vec{r},t)}{\partial t}=-k_{AHL}\cdot C_{AHL}(\vec{r},t)
 
                        \;\;\; \text{(3a)} $$
 
                        $$ \frac{\partial C_{leucine}(\vec{r},t)}{\partial t}=-k_{leucine}\cdot C_{leucine}(\vec{r},t)
 
                        \;\;\; \text{(3b)} $$
 
</p>
 
</div>
 
<p>
 
                        Putting it all together, we obtain (4), both for AHL and leucine.
 
</p>
 
<div id="eq4">
 
<p>
 
                        $$ \frac{\partial C(\vec{r},t)}{\partial t}=D \cdot \nabla^2 C(\vec{r},t)+\alpha \cdot \rho_A(\vec{r},t)-k\cdot
 
                        C(\vec{r},t) \;\;\; \text{(4)} $$
 
</p>
 
</div>
 
<p>
 
                        Note that these equations have exactly the same form as the equations for AHL and leucine
 
                        in the colony level model. The crucial difference however lies in the
 
                        calculation of the density of cells of type A. In contrast to the colony level
 
model, in this model the cell density is
 
not calculated explicitly with a PDE and is therefore
 
                        not trivially known. Therefore a method to extract a density field from a
 
                        spatial distribution of agents is necessary. This is addressed in the
 
                        subparagraph below on coupling.
 
                    </p>
 
 
</div><!-- part -->
 
 
<div class="center">
 
<div class="togglebar">
 
<div class="toggleone">
 
                        <h2>
 
                            Box 1: Heat Equation
 
                        </h2>
 
</div>
 
<div id="toggleone">
 
<div class="widebox">
 
                  <h2> Heat Equation </h2>
 
                  <p>
 
  The left-hand side of the heat equation (1) represents the rate of
 
  accumulation of a chemical, while the right-hand side is proportional
 
  to the Laplacian
 
  of its concentration field, which is a second-order differential operator.
 
  This equation can be easily understood by
 
  considering a one-dimensional concentration profile, as shown in Figure 1: if the concentration
 
  can be approximated as a convex parabolic function, the second derivative
 
  is positive and therefore the rate of accumulation is positive (i.e. more
 
  accumulation). If on the other hand the concentration resembles a concave
 
  parabolic function, the second derivative is negative and the rate of
 
  accumulation as well (i.e. depletion). A special case occurs when the
 
  concentration profile takes on a linear form. Everything that moves into
 
  the point goes out at the other side and as a result there is no net accumulation
 
  over time.
 
                  </p>
 
 
                  <div class="center">
 
                        <div id="image1">
 
                        <a class="example-image-link" href="https://static.igem.org/mediawiki/2015/0/0a/KU_Leuven_heatDiagram.png" data-lightbox="example-set" data-title="Illustration of the heat equation"><img class="example-image" src="https://static.igem.org/mediawiki/2015/0/0a/KU_Leuven_heatDiagram.png" alt="Illustration of the heat equation" width="45%" height="45%"></a>
 
<h4><div id=figure1>Figure 1</div> Illustration of the heat equation. Click to enlarge.</h4>
 
</div><!-- image1 -->
 
</div><!-- center -->
 
 
<p>
 
In the videoplayer below we demonstrate how the heat equation smoothes out
 
an initially heterogeneous concentration profile. When only diffusion is
 
acting on the system, it will always evolve to a uniformly flat
 
concentration profile, regardless of the initial conditions.
 
</p>
 
<div class="center">
 
                    <!-- third Videobox start-->
 
                    <video id="video3" preload="auto" width="50%" tabindex="0" controls="" type="video/mp4">
 
                      <source type="video/mp4" src="https://static.igem.org/mediawiki/2015/6/6e/KU_Leuven_EpanechnikovHeat.mp4">
 
                      Sorry, your browser does not support HTML5 video.
 
                    </video>
 
                    <br/>
 
            <button type="button" onclick="Set8()">Epanechnikov</button>
 
            <button type="button" onclick="Set9()">iGEM</button>
 
            <script>
 
            function Set8() {
 
            document.querySelector("#video3 > source").src = "https://static.igem.org/mediawiki/2015/6/6e/KU_Leuven_EpanechnikovHeat.mp4"
 
            document.querySelector("#video3").load();
 
            document.querySelector("#video3").play();
 
            }
 
            function Set9() {
 
            document.querySelector("#video3 > source").src = "https://static.igem.org/mediawiki/2015/9/9f/KU_Leuven_iGEMHeat.mp4"
 
            document.querySelector("#video3").load();
 
            document.querySelector("#video3").play();
 
            }     
 
            </script>
 
</div><!-- center -->
 
            <!-- video end -->
 
<div class="whiterow"></div>
 
</div><!-- widebox -->
 
</div><!-- toggleone -->
 
</div><!-- togglebar -->
 
 
</div><!-- center, end of toggles -->
 
 
                <div class="part">
 
<h2>Agents</h2>
 
                    <p>
 
                        To model bacteria movement on the other hand, we used an agent-based model that
 
                        explicitly stored individual bacteria as agents. Spatial coordinates are
 
                        associated with each agent, specifying their location. After solving the
 
                        equation of motion for all agents based on their environment, these coordinates
 
                        are updated at every timestep. In principle, Newton’s second law of motion has
 
                        to be solved for all bacteria. However, since bacteria live in a low Reynolds
 
                        (high friction) environment, the inertia of the bacteria can be neglected. This
 
                        is because an applied force will immediately be balanced out by an opposing
 
                        frictional force, with no noticeable acceleration or deceleration phase taking
 
                        place.                                       
 
                        This eliminates the inertial term and simplifies Newton’s second law to
 
                        (5).
 
</p>
 
<div id="eq5">
 
<p>
 
                        $$  \frac{d^2 \vec{r}(t)}{dt^2}=\sum_{i} \vec{F}_{applied,i}-\gamma \cdot
 
                        \frac{d \vec{r}(t)}{dt}=0 $$
 
                        $$\Rightarrow \frac{d \vec{r}(t)}{dt}=\frac{1}{\gamma}
 
                        \cdot \sum_{i} \vec{F}_{applied,i} \;\;\; \text{(5)} $$
 
</p>
 
</div>
 
<p>
 
                        Basically, the velocity can be calculated as the sum of all applied
 
forces times divided by a frictional coefficient.
 
For more info about the Reynolds number and
 
“life at low Reynolds number”, we refer to box 2. In the following
 
paragraphs we will investigate the different forces acting on the bacteria
 
and ultimately superimpose them to obtain the final equation of motion.
 
                    </p>
 
</div><!-- part -->
 
 
<div class="center">
 
<div class="togglebar">
 
<div class="toggle14">
 
                        <h2>
 
                            Box 2: Life at Low Reynolds Number
 
                        </h2>
 
</div>
 
<div id="toggle14">
 
<div class="widebox">
 
                  <h2> Life at Low Reynolds Number </h2>
 
                  <p>
 
  The Reynolds number in fluid mechanics is a dimensionless number that characterizes
 
  different flow regimes. It is most commonly used to determine whether laminar
 
  or turbulent flow will take place in a hydrodynamic system (Figure 2).
 
  </p>
 
 
                  <div class="center">
 
<div id="image2">
 
<a class="example-image-link"
 
href="https://static.igem.org/mediawiki/2015/6/61/KU_Leuven_laminar_vs_turbulent_flow.png"
 
data-lightbox="example-set"
 
data-title="Flow regimes for different Reynolds numbers">
 
<img class="example-image"
 
src="https://static.igem.org/mediawiki/2015/6/61/KU_Leuven_laminar_vs_turbulent_flow.png"
 
alt="Flow regimes for different Reynolds numbers"
 
width="70%"
 
></a>
 
<h4><div id=figure2>Figure 2</div>
 
Flow regimes for different Reynolds numbers. Click to enlarge.</h4>
 
</div><!-- image2 -->
 
</div><!-- center -->
 
 
  <p>
 
  In general however, it quantifies the ratio of inertial forces to viscous forces
 
  and is defined as (B2.1).
 
  </p>
 
 
<div id="eq20">
 
<p>
 
$$
 
Re=\frac{\text{inertial forces}}{\text{viscous forces}}
 
=\frac{\rho v^2 L^2}{\eta v L}
 
=\frac{\rho v L}{\eta} \;\;\; \text{(B2.1)}
 
$$
 
</p>
 
</div>
 
 
<p>
 
For example, if the Reynolds number is high, the inertia of fluids in motion dominate and
 
turbulent flow will occur. When it is low however, the viscous forces dampen
 
the kinetic energy of fluid particles and stabilize the flow profile, ultimately
 
achieving a regular, laminar flow.
 
<br/><br/>
 
The reason we mention the Reynolds number however is not to study the flow of fluids, but
 
to characterize the behavior of swimming objects inside of a stationary fluid.
 
When applied to objects in a fluid, the Reynolds number tells us whether their inertia
 
can be neglected or not. For bacteria, the characteristic length $L$ is on the order
 
of micrometers, which is quite small. Therefore the Reynolds number is always
 
very small and hence viscous forces dominate the motion of bacteria. This is what allows
 
us to eliminate the inertial term in Newton's second law and greatly simplify the
 
equation of motion.
 
<br/><br/>
 
To justify this simplification we will go through a numerical example.
 
Take a bacterium of size $L = 2 a = 1 \cdot \mu m$, swimming through water with a speed of
 
$v=20 \cdot \mu m/s$. Water has a density of $\rho = 1 \; 000 \cdot kg/m^3$
 
and a viscosity of $\eta = 0.001 \cdot Pa/s$.
 
</p>
 
 
                  <div class="center">
 
<div id="image3">
 
<a class="example-image-link"
 
href="https://static.igem.org/mediawiki/2015/5/5c/KU_Leuven_diagram_bacterium_reynolds.png"
 
data-lightbox="example-set"
 
data-title="Sketch of swimming bacterium">
 
<img class="example-image"
 
src="https://static.igem.org/mediawiki/2015/5/5c/KU_Leuven_diagram_bacterium_reynolds.png"
 
alt="Sketch of swimming bacterium"
 
width="70%"
 
></a>
 
<h4><div id=figure2>Figure 3</div>
 
Sketch of swimming bacterium. Click to enlarge.</h4>
 
</div><!-- image3 -->
 
</div><!-- center -->
 
 
<p>
 
Putting these values in (B2.1) yields
 
a Reynolds number of $Re = 2 \cdot 10^{-5} << 1$. Clearly we are in an extremely
 
low Reynolds number regime. To show what this really means, suppose the bacterium
 
stops propelling itself. How long will it continue to move relying only on its inertia?
 
Assuming Stokes' law, we obtain (B2.2) as the equation of motion.
 
</p>
 
 
 
<div id="eq21">
 
<p>
 
$$
 
F=m\cdot \frac{d^2x}{dt^2} \;\;\; \text{(B2.2a)}
 
$$
 
$$
 
-6 \pi \eta a  v=m\cdot \frac{dv}{dt} \;\;\; \text{(B2.2b)}
 
$$
 
</p>
 
</div>
 
 
<p>
 
Solving this differential equation yields (B2.3).
 
</p>
 
 
<div id="eq22">
 
<p>
 
$$
 
v=v_0 \cdot \text{exp} \Bigg[
 
-\frac{6 \pi \eta a}{m} \cdot (t-t_0)
 
\Bigg]
 
\;\;\; \text{(B2.3)}
 
$$
 
</p>
 
</div>
 
 
<p>
 
The characteristic time constant is $\tau = 6\pi\eta a/m \approx 0.1 \cdot \mu s$,
 
from which we calculate that the total distance traversed is $\Delta x < 2 \cdot pm$.
 
This coasting distance is 6 orders of magnitude smaller than its size. Moreover,
 
the time spent coasting is extremely short. Thus, once the
 
bacterium stops propelling itself, it is safe to assume that whatever kinetic energy
 
it had is immediately absorbed by hydrodynamic friction, instantly halting the bacterium.
 
Therefore, we can neglect the inertial term in Newton's second law (5).
 
</p>
 
 
<div class="whiterow"></div>
 
</div><!-- widebox -->
 
</div><!-- toggleone -->
 
</div><!-- togglebar -->
 
 
              <div class="togglebar">
 
                <div class="toggletwo">
 
                    <h2>
 
                        Stochastic Differential Equation
 
                    </h2>
 
</div><!-- toggletwo -->
 
              <div id="toggletwo">
 
  <p>
 
  When we try to calculate the physical “chemotactic force”, propelling bacteria towards
 
  chemoattractants or away from chemorepellents, we find that it is not easily
 
  measured nor derived. Therefore, as a workaround
 
  we base the equation of motion
 
                        on (6), a stochastic differential equation (SDE) that describes the motion
 
                        of a single particle in a N-particle system that is governed by a Keller-Segel
 
type PDE in the limit of $N \rightarrow \infty$. This Keller-Segel type PDE (7)
 
describes the evolution of a bacteria density field $n$ moving
 
towards some chemoattractant $S$.
 
</p>
 
<div id="eq6">
 
<p>
 
$$ d\vec{r}_i(t)=\chi (S)
 
\cdot \nabla S(\vec{r},t)\cdot dt + \sqrt{2 \cdot
 
\mu}\cdot d\vec{W} \;\;\; \text{(6a)} $$
 
$$ d\vec{r}_i(t)=\mu \cdot \frac{\kappa}{S(\vec{r},t)}
 
\cdot \nabla S(\vec{r},t)\cdot dt + \sqrt{2 \cdot
 
\mu}\cdot d\vec{W} \;\;\; \text{(6b)} $$
 
</p>
 
</div>
 
<p>
 
</p>
 
<div id="eq7">
 
<p>
 
$$ \frac{\partial n(\vec{r},t)}{\partial t}=\mu \cdot \nabla^2 n
 
- \nabla (n \cdot \chi(S)
 
\cdot\nabla S(\vec{r},t)) \;\;\; \text{(7a)} $$
 
 
$$ \frac{\partial n(\vec{r},t)}{\partial t}=\mu \cdot \nabla^2 n
 
- \nabla (n \cdot \mu \cdot \frac{\kappa}{S(\vec{r},t)}
 
\cdot\nabla S(\vec{r},t)) \;\;\; \text{(7b)} $$
 
</p>
 
</div>
 
<p>
 
Put differently, when infinitely many particles obey (6), they will exhibit
 
Keller-Segel type spatial dynamics as in (7). In a sense,
 
                        we’re using a “reverse-engineered” particle equation that corresponds to the
 
                        Keller-Segel field equation. A detailed theoretical treatment of (6) is
 
                        outside the scope of this model description because it contains a stochastic
 
                        variable. The traditional rules of calculus do not apply anymore for stochastic
 
                        differential equations and a different mathematical theory called Itô calculus
 
is required. It is sufficient to say that the second term containing
 
$dW$ accounts
 
                        for Brownian motion in the form of random noise added to the displacement of the
 
                        agents, causing them to diffuse, and that it is governed by the diffusion
 
coefficient $\mu$.
 
<br/>
 
<br/>
 
  The first term in (6) on the other hand is easily understood
 
                        as an advective or drift term (net motion) dependent on S, pushing the agents
 
                        along a positive gradient (for negative chemotaxis the sign is reversed). The
 
                        chemotactic drift hence points towards an increasing concentration of the
 
                        chemoattractant. The advective properties are governed by the chemotactic
 
sensitivity function $\chi (S)$. For our model we used a function of the
 
form (8).
 
</p>
 
<div id="eq8">
 
<p>
 
$$ \chi(S)=\mu \cdot \frac{\kappa}{S(\vec{r},t)} \;\;\; \text{(8)} $$
 
</p>
 
</div>
 
<p>
 
                        The first important thing to note is that we assume
 
                        $\chi (S)$ to be proportional to $1/S$. This is because Keller and Segel proved that
 
                        their corresponding PDE model only yields travelling wave solutions if $\chi (S)$
 
contains such a singularity at some critical concentration $S_{crit}$,
 
and multiplying by
 
                        $1/S$ is the simplest way to introduce a singularity at $S_{crit} = 0$. Secondly, the
 
                        proportionality constant is composed of two factors, namely the bacterial
 
diffusion coefficient $\mu$ and chemotactic sensitivity constant $\kappa$.
 
This is done for
 
                        two reasons. Firstly, when $\mu$ is lowered, both chemotactic and random motion is
 
                        reduced, which emulates the state of inactivated motility due to high or low
 
                        concentrations of AHL. Secondly, defining a separate chemotactic sensitivity
 
                        constant allows us to examine the effect of the relative strength of chemotaxis
 
versus random motion on pattern formation. Conceptually, this term can be
 
viewed as a sort
 
of chemotactic force, defined with regards to a mobility coefficient $\mu$
 
instead of a frictional coeffient, which are the inverse of each other.
 
                    </p>
 
</div><!-- toggletwo -->
 
</div><!-- togglebar -->
 
 
        <div class="togglebar">
 
            <div class="togglethree">
 
                    <h2>
 
                        Cell-Cell Interactions
 
                    </h2>
 
</div><!-- togglethree -->
 
            <div id="togglethree">
 
                    <p>
 
                        In addition to chemotaxis and diffusion, cell-cell interactions play an
 
                        important role in pattern formation and also need to be modeled. Bacteria have
 
                        finite size and therefore multiple bacteria cannot occupy the same space.
 
                        Moreover, an important mechanism in our system is the aggregation of cells A due
 
                        to the sticky adhesin protein membrane. To take these mechanisms into account we
 
                        modeled two types of cell-cell interactions: the purely repulsive interaction of
 
                        cell B with another cell B and with cell A, and the attractive-repulsive
 
                        interaction of cell A with another cell A. The interaction between two cells is
 
                        usually expressed by a potential energy curve defined over the distance between
 
the centers of mass of the two cells (Figure 2).
 
</p>
 
<div class="center">
 
<div id="image2">
 
                        <a class="example-image-link" href="https://static.igem.org/mediawiki/2015/e/e6/KU_Leuven_CellCellInt.png" data-lightbox="example-set" data-title="Illustration of cell-cell interaction potential"><img class="example-image" src="https://static.igem.org/mediawiki/2015/e/e6/KU_Leuven_CellCellInt.png" alt="Illustration of cell-cell interaction potential" width="45%" height="45%"></a>
 
                        <h4><div id=figure2>Figure 2</div>Cell-cell interaction potential. Click to enlarge.</h4>
 
</div><!-- image2 -->
 
</div><!-- center -->
 
<p>
 
Note that the potential energy
 
                        remains constant after a certain distance, which means that the cells stop
 
                        interacting. Also, as two cells move closer together, they hit a wall where the
 
                        potential energy curve abruptly goes to infinity. The reason for this is that
 
                        two cells cannot occupy the same space and therefore smaller intercellular
 
                        distances are not allowed. Implementing this theoretical potential is however
 
not possible because the displacement of bacteria is a stochastic process
 
and the bacteria could randomly jump beyond
 
                        the potential wall, where the force is ill defined. Therefore, we’ve decided
 
                        to define a piecewise quadratic potential (9),
 
                        which results in a piecewise
 
                        linear force that resembles Hooke’s law, but with three different “spring
 
constants” acting in different intervals of intercellular distances (10),
 
as illustrated in Figure 3. The force is defined with respect to
 
the unit vector pointing towards the other cell, meaning that a positive
 
force corresponds to an attractive force and vice-versa.
 
</p>
 
 
<div class="center">
 
<div id="image3">
 
                        <a class="example-image-link" href="https://static.igem.org/mediawiki/2015/8/89/KU_Leuven_potential_force_curve.png" data-lightbox="example-set" data-title="Cell-cell interaction potential and force curves"><img class="example-image" src="https://static.igem.org/mediawiki/2015/8/89/KU_Leuven_potential_force_curve.png" alt="Cell-cell interaction potential and force curves" width="45%" height="45%"></a>
 
                        <h4><div id=figure3>Figure 3</div>Cell-cell interaction potential and force curves. Click to enlarge.</h4>
 
</div><!-- image2 -->
 
</div><!-- center -->
 
<p>
 
 
<div id="eq9">
 
<p>
 
 
                        $$ E_{p,attr-rep}(r_{ij})=\left\{\begin{matrix} 0
 
& R_{cutoff}\leq r_{ij}\\
 
  \frac{1}{2}\cdot k_3 \cdot(r_{ij}-R_0)^2+C_1
 
  & R_0 \leq r_{ij} < R_{cutoff} \\
 
  \frac{1}{2}\cdot k_2 \cdot(r_{ij}-R_0)^2+C_1
 
  & \frac{R_0}{2} \leq r_{ij} < R_0\\
 
  \frac{1}{2}\cdot k_1 \cdot(r_{ij}-\frac{k_1+k_2}{k_1}\cdot \frac{R_0}{2})^2+C_2
 
  & 0 \leq r_{ij} < \frac{R_0}{2}
 
  \end{matrix}\right. \;\;\; \text{(9a)}$$
 
 
  $$ E_{p,rep}(r_{ij})=\left\{\begin{matrix} 0
 
  & R_0\leq r_{ij}\\
 
  \frac{1}{2}\cdot k_2 \cdot(r_{ij}-R_0)^2
 
  & \frac{R_0}{2} \leq r_{ij} < R_0\\
 
  \frac{1}{2}\cdot k_1 \cdot(r_{ij}-\frac{k_1+k_2}{k_1}\cdot \frac{R_0}{2})^2
 
  & 0 \leq r_{ij} < \frac{R_0}{2}
 
  \end{matrix}\right. \;\;\; \text{(9b)}$$
 
 
</p>
 
</div>
 
 
<div id="eq10">
 
<p>
 
$$ \vec{F}_{i,attr-rep}(r_{ij})=
 
\frac{\partial E_{p,attr-rep}(r_{ij})}{\partial r_{ij}} \cdot \vec{e}_{ij}=
 
\left\{\begin{matrix}
 
\vec{0} & R_{cutoff}\leq r_{ij}\\
 
k_3 \cdot(r_{ij}-R_0) \cdot \vec{e}_{ij}
 
& R_0 \leq r_{ij} < R_{cutoff} \\
 
k_2 \cdot(r_{ij}-R_0) \cdot \vec{e}_{ij}
 
&\frac{R_0}{2} \leq r_{ij} < R_0\\
 
k_1 \cdot(r_{ij}-\frac{k_1+k_2}{k_1}\cdot \frac{R_0}{2}) \cdot \vec{e}_{ij}
 
&  0 \leq r_{ij} < \frac{R_0}{2}
 
\end{matrix}\right. \;\;\; \text{(10a)}
 
                        $$
 
 
$$ \vec{F}_{i,rep}(r_{ij})=
 
\frac{\partial E_{p,rep}(r_{ij})}{\partial r_{ij}} \cdot \vec{e}_{ij}=
 
\left\{\begin{matrix}
 
                          \vec{0} & R_0\leq r_{ij}\\
 
  k_2 \cdot(r_{ij}-R_0) \cdot \vec{e}_{ij}
 
  & \frac{R_0}{2} \leq r_{ij} < R_0\\
 
  k_1 \cdot(r_{ij}-\frac{k_1+k_2}{k_1}\cdot \frac{R_0}{2}) \cdot \vec{e}_{ij}
 
  &  0 \leq r_{ij} < \frac{R_0}{2}
 
  \end{matrix}\right. \;\;\; \text{(10b)}
 
                        $$
 
</p>
 
</div>
 
<p>
 
 
Between A type cells, there is a region of attraction
 
$(R_0 \leq r_{ij} < R_{cutoff})$,
 
where the force points towards the other cell,
 
hence moving them closer together.
 
In the repulsive domain $(r_{ij} < R_0)$,
 
two regions were defined, emulating lower
 
repulsive forces $(\frac{R_0}{2} \leq r_{ij} < R_0)$
 
and higher repulsive forces due to a higher spring
 
constant when the cells are even closer $(r_{ij} < \frac{R_0}{2})$.
 
For the purely repulsive interaction
 
scheme there is no attraction and therefore
 
the spring constant for $R_0 \leq r_{ij}$ is zero.
 
More details about the implementation of the
 
cell-cell interaction scheme, more specifically
 
regarding the nearest-neighbor search algorithm,
 
can be found in the paragraph on the agent-based
 
module below.
 
                    </p>
 
</div><!-- togglethree -->
 
</div><!-- togglebar -->
 
 
    <div class="togglebar">
 
        <div class="togglefour">
 
<h2>Equation of Motion</h2>
 
</div>
 
<div id="togglefour">
 
                    <p>
 
                        Now we are ready to construct the equation of motion for cell type A and B as a
 
superposition of an adapted Keller-Segel SDE (6)
 
and cell-cell interaction forces (10),
 
                        yielding (11).
 
</p>
 
<div id="eq11">
 
<p>
 
$$  d\vec{r}_{A_i}(t)= \sqrt{2 \cdot \mu_A (H)}\cdot d\vec{W}
 
+ \frac{1}{\gamma}\cdot\Bigg[ \sum^{A
 
\backslash \{ A_i\}}_j \frac{dE_{p,attr-rep}(r_{ij})}{dr_{ij}}
 
\cdot \vec{e}_{ij}
 
+\sum^{B}_j\frac{dE_{p,rep}(r_{ij})}{dr_{ij}}\cdot
 
\vec{e}_{ij} \Bigg]\cdot dt \;\;\; \text{(11a)}
 
                        $$
 
                        $$
 
d\vec{r}_{B_i}(t)=
 
\sqrt{2 \cdot \mu_B(H)}\cdot d\vec{W} +
 
\chi(L,H) \cdot \nabla L(\vec{r},t)\cdot dt +
 
\frac{1}{\gamma}\cdot\Bigg[ \sum^{A\cup B\backslash
 
\{ B_i\}}_j \frac{dE_{p,rep}(r_{ij})}{dr_{ij}}\cdot \vec{e}_{ij} \Bigg]\cdot dt
 
\;\;\; \text{(11b)}
 
                        $$
 
</p>
 
</div>
 
 
<div id="eq12">
 
<p>
 
                        $$
 
\chi(L,H)= -\mu_{B}(H)
 
\cdot \frac{\kappa}{L(\vec{r},t)} \;\;\; \text{(12)}
 
                        $$
 
</p>
 
</div>
 
 
<div id="eq13">
 
<p>
 
                        $$
 
\mu_A(H)=\left\{\begin{matrix}
 
\mu_{A,high} & H(\vec{r},t) < H_{A,threshold}\\
 
\mu_{A,low} & H(\vec{r},t) \geq H_{A,threshold}
 
\end{matrix}\right.
 
\;\;\; \text{(13a)}
 
                        $$
 
                        $$
 
\mu_B(H)=\left\{\begin{matrix}
 
\mu_{B,high} & H(\vec{r},t) < H_{B,threshold}\\
 
\mu_{B,low} & H(\vec{r},t) \geq H_{B,threshold}
 
\end{matrix}\right.
 
\;\;\; \text{(13b)}
 
                        $$
 
</p>
 
</div>
 
 
<p>
 
                        Bacteria of type A are not attracted nor repelled by leucine,
 
                        so the chemotactic term falls away. All cell-cell forces are summed up to find a
 
                        net force, taking into account the two different potentials due to the different
 
                        interaction types. As discussed before, this net force times a constant yields
 
                        the velocity due to that force, which is then multiplied by $dt$ to obtain the
 
                        displacement. For B type cells, the chemotactic term models the repulsive
 
                        chemotaxis away from leucine. The chemotactic sensitivity function (12) has a
 
                        negative sign signifying that B type cells are repelled by leucine. The cell-cell
 
                        interaction term in this case is simpler because B type cells only interact
 
repulsively.
 
<br/>
 
<br/>
 
Note that the diffusion coefficient of cell types A and B (13) switches
 
                        based on the local concentration of AHL relative to a threshold AHL value, which
 
simulates the dependency of cellular motility on AHL.
 
For B type cells the cellular motility depends explicitly on AHL due to the
 
synthetic genetic circuit we have built into them. On the other hand,
 
in our model the motility of A type cells should not depend on AHL directly,
 
but since high concentrations of AHL are caused by high densities of
 
aggregated A type cells, the bacteria typically will not be motile
 
in high AHL
 
concentrations because they stick to neighboring cells.
 
<br/>
 
<br/>
 
The agent-based module is
 
                        now fully defined but one crucial issue was skipped: AHL and leucine
 
                        concentrations are calculated using PDEs and are therefore only known at grid
 
                        points. Agents on the other hand reside in the space between grid points and
 
                        require local concentrations as inputs to calculate their next step. This
 
                        problem is part of the coupling aspect in our hybrid modeling framework and is
 
                        discussed below.
 
                    </p>
 
         
 
</div><!-- togglefour -->
 
</div><!-- togglebar -->
 
</div><!-- center -->
 
 
<div class="part">
 
<h2>Coupling</h2>
 
<p>
 
At this point, both the PDE module and agent-based module have been established,
 
but the issue remains that the individual modules take inputs and return outputs which are not
 
directly compatible. In the framework of PDEs, entities are described by
 
density/concentration fields.
 
In the agent-based module however, entities are represented by discrete
 
objects with exact spatial locations. The key challenge posed by the hybrid model
 
is to integrate these different approaches by bidirectionally converting and
 
exchanging information between the modules. Coupling the modules is the essential component
 
that makes the hybrid model hybrid. It allows us to leverage the advantages
 
of both modeling approaches, while circumventing their drawbacks. But although
 
hybrid modeling opens up many new avenues for novel modeling methods,
 
it comes with its own diverse set of issues and peculiarities that need to be
 
addressed before it can be successfully applied.
 
<br/><br/>
 
In the following paragraphs the basic scheme for coupling
 
the PDE and agent-based modules in our model is introduced, after which the theoretical
 
treatment of our hybrid model is complete. The difficulties that arise from linking the
 
modules together are discussed further below in the subsection on matching.
 
 
</p>
 
</div><!-- part -->
 
 
<div class="center"><!-- start of toggles -->
 
<div class="togglebar">
 
<div class="togglefive">
 
<h2>
 
Agent-Based to PDE
 
</h2>
 
</div><!-- togglefive -->
 
<div id="togglefive">
 
                    <p>
 
                        As described above, the agents’ effect in the PDE is modeled as a source term
 
                        that is proportional to the agent density. This approach is essentially the same
 
                        approach taken in the colony level model for the bacterial production of AHL and
 
                        leucine. However, in the colony level model the bacteria density is explicitly
 
                        calculated at the grid points, while the agent-based model essentially considers
 
a set of points scattered in space.
 
A simple first-order approach would be to determine
 
                        the closest grid point to any agent and simply increment a counter belonging to
 
                        that grid point. This results in a histogram, which can be used directly to
 
                        represent the agent density. However, the resulting density is a blocky,
 
                        nonsmooth function which poorly represents the underlying agent distribution.
 
                        The effect of a single agent is artificially confined to a single grid point,
 
                        while in reality an agent’s influence could reach much further than a single
 
                        grid point. The shape of a histogram is also very dependent on the bin size,
 
                        which in this case corresponds to the grid spacing so it cannot be independently
 
                        tuned. To decouple grid spacing and agent density and achieve a smoother density
 
                        function, we made use of a more sophisticated technique called kernel density
 
                        estimation (KDE).
 
                    </p>
 
                    <p>
 
                        KDE is used in statistics to estimate the probability density of a set discrete
 
                        data derived from a random process. The basic idea consists of defining a
 
                        kernel function that represents the density of a single data point, then
 
                        centering kernel functions on every data point and summing them all up to
 
                        achieve a smooth overall density function, as demonstrated in the Figure 4. </p>
 
 
<div class="center">
 
<div id="image4">
 
                        <a class="example-image-link" href="https://static.igem.org/mediawiki/2015/9/9c/KU_Leuven_KernelSum.png"
 
data-lightbox="example-set" data-title="Kernel density estimation">
 
<img class="example-image"
 
                            src="https://static.igem.org/mediawiki/2015/9/9c/KU_Leuven_KernelSum.png"
 
                            alt="Kernel sum" width="50%" ></a>
 
                        <h4><div id=figure4>Figure 4</div> Kernel density estimation. Click to enlarge.</h4>
 
</div><!-- image1 -->
 
</div><!-- center -->
 
 
                      <p>
 
                      This
 
                        kernel function can be anything as long is it continuous, symmetric and
 
                        integrates to 1, since it represents one data point or one agent in our case.
 
                        Some of the most common kernel functions include Gaussian kernels, triangular
 
                        kernels and Epanechnikov kernels (14), which are shown in Figure 5.
 
</p>
 
 
                      <div class="center">
 
                        <div id="image5">
 
<a class="example-image-link"
 
href="https://static.igem.org/mediawiki/2015/8/88/KU_Leuven_variousKernelTypes.png"
 
data-lightbox="example-set"
 
data-title="Gaussian, triangular and Epanechnikov kernel functions">
 
<img class="example-image"
 
src="https://static.igem.org/mediawiki/2015/8/88/KU_Leuven_variousKernelTypes.png"
 
alt="Gaussian, triangular and Epanechnikov kernel functions"
 
width="60%"></a>
 
<h4><div id=figure1>Figure 5</div>
 
Gaussian, triangular and Epanechnikov kernel functions.
 
Click to enlarge.</h4>
 
</div><!-- image5 -->
 
</div><!-- center -->
 
 
<div id="eq14">
 
<p>
 
$$
 
K_{Gaus}(x)=\frac{1}{\sqrt{2 \cdot \pi}}
 
\cdot e^{-\frac{1}{2}x^2} \;\;\; \text{(14a)}
 
$$
 
$$
 
K_{tri}(x)=\left\{\begin{matrix}
 
0 & 1 < |x| \\
 
1-|x| & |x| \leq 1
 
\end{matrix}\right. \;\;\; \text{(14b)}
 
$$
 
$$
 
K_{Epa}(x)=\left\{\begin{matrix}
 
0 & 1 < |x| \\
 
\frac{3}{4} \cdot (1-x^2) & |x| \leq 1
 
\end{matrix}\right. \;\;\; \text{(14c)}
 
$$
 
</p>
 
</div>
 
<p>
 
 
In practice, scaled versions of standard kernel functions are used,
 
which are of the form (15).
 
 
</p>
 
<div id="eq15">
 
<p>
 
$$
 
K_h(x)=\frac{1}{h} \cdot K \bigg( \frac{x}{h} \bigg) \;\;\; \text{(15)}
 
$$
 
</p>
 
</div>
 
<p>
 
                        Importantly, the scaled functions inherit the kernel function
 
                        properties, but are either broader or narrower. The degree to which the shape of
 
                        a kernel function is stretched or squeezed depends on the scaling factor h,
 
which is why it is called the <i>bandwidth</i>, see Figure 6.
 
</p>
 
                      <div class="center">
 
                        <div id="image6">
 
                        <a class="example-image-link" href="https://static.igem.org/mediawiki/2015/5/52/KU_Leuven_EpanechnikovDifferentBandwidth.png" data-lightbox="example-set" data-title="Epanechnikov kernel using different bandwidths"><img class="example-image" src="https://static.igem.org/mediawiki/2015/5/52/KU_Leuven_EpanechnikovDifferentBandwidth.png" alt="Epanechnikov kernel using different bandwidths" width="45%"></a>
 
                        <h4><div id=figure6>Figure 6</div>Epanechnikov kernels using different bandwidths. Click to enlarge.</h4>
 
</div><!-- image3 -->
 
</div><!-- center -->
 
<p>
 
This parameter gives us the
 
                        freedom to define how far the influence of an agent reaches and how smooth the
 
resulting density function looks like. An example is given in Figure 7, where
 
100 agents were sampled from a normal distribution. The red dotted line
 
indicates the true underlying distribution, while the blue line is the
 
kernel density estimation calculated using different bandwidths. The left graph
 
shows a undersmoothed (small bandwidth) KDE, which oscillates wildly.
 
The right graph on the other hand is oversmoothed (large bandwidth),
 
leading to a misleadingly wide KDE.
 
</p>
 
                      <div class="center">
 
                        <div id="image7">
 
<a class="example-image-link" href="https://static.igem.org/mediawiki/2015/1/15/KU_Leuven_bandwidth_density.png" data-lightbox="example-set" data-title="Kernel density estimation using different bandwidths"><img class="example-image" src="https://static.igem.org/mediawiki/2015/1/15/KU_Leuven_bandwidth_density.png" alt="Kernel density estimation using different bandwidths" width="70%"></a>
 
 
                        <h4><div id=figure7>Figure 7</div>Kernel density estimation using different bandwidths. Click to enlarge.</h4>
 
</div><!-- image3 -->
 
</div><!-- center -->
 
<p>
 
In summary, using kernel density estimation we can express the agent density
 
in the form of (16), providing the appropriate input for the PDE model.
 
</p>
 
<div id="eq16">
 
<p>
 
$$
 
\rho(x)=\beta \cdot \sum^N_{i} K_h(x-x_i)
 
\;\;\; \text{(16)}
 
$$
 
</p>
 
</div>
 
<p>
 
The extension to higher dimensions is called multivariate kernel density
 
estimation and is rather non-trivial since the bandwidth is then defined
 
as a matrix instead of a scalar, which allows for many more variations
 
on the same basis kernel function. A detailed analysis of multivariate
 
kernel density estimation will not be given here. It suffices to
 
say that for our 2-D hybrid model we used the simplest
 
possible bandwidth matrix, which is the identity matrix times a constant.
 
</p>
 
 
</div><!-- togglefive -->
 
</div><!-- togglebar -->
 
 
<div class="togglebar">
 
<div class="togglesix">
 
<h2>
 
PDE to Agent-Based
 
</h2>
 
</div>
 
<div id="togglesix">
 
                    <p>
 
The final component of our hybrid model is the mapping
 
of the PDE model to the agent-based model.
 
The latter model works with discrete objects that
 
have continuous coordinates, meaning that they can be located
 
                  at any point of the domain.
 
  As we have seen in (11), the agents need the local concentration
 
  of AHL and leucine, as well as the gradient of leucine in
 
  order to update their positions. In the PDE model however,
 
  the domain is discretized into a grid and concentrations
 
  are only defined at grid points. Therefore, in order to
 
  transfer information from the PDE model to the agent-based
 
  model we need to translate these grid values into values
 
  for any given position within the domain. We achieved this
 
  for the 1-D hybrid model by taking the two nearest grid points
 
  to any agent and employing linear interpolation to
 
  derive an approximate local field value. Similarly, for the 2-D
 
  hybrid model we took the four nearest grid points
 
                  and employed bilinear interpolation, as explained in box 3.
 
  With the end result (B3.3), the concentration or gradient of AHL
 
  or leucine can be derived at any point of the domain. Therefore,
 
  the agent-based module has access to the information that it needs
 
  from the PDE module.
 
                    </p>
 
</div><!-- togglesix -->
 
</div><!-- togglebar -->
 
 
<div class="togglebar">
 
<div class="toggleseven">
 
<h2>
 
Box 3: Bilinear Interpolation
 
</h2>
 
</div>
 
<div id="toggleseven">
 
<div class="widebox">
 
<h2> Bilinear Interpolation </h2>
 
<p>
 
Interpolation is used to obtain a continuous function
 
when a finite amount of data points are known. Piecewise
 
linear interpolation is the most simple scheme to interpolate
 
1-dimensional data.
 
</p>
 
 
  <div class="center">
 
<div id="image8">
 
<a class="example-image-link"
 
href="https://static.igem.org/mediawiki/2015/9/9f/KU_Leuven_linear_interpolation.png"
 
data-lightbox="example-set"
 
data-title="Illustration of linear interpolation">
 
<img class="example-image"
 
src="https://static.igem.org/mediawiki/2015/9/9f/KU_Leuven_linear_interpolation.png"
 
alt="Illustration of linear interpolation"
 
width="60%"
 
></a>
 
<h4>
 
<div id=figure8>Figure 8</div>
 
Illustration of linear interpolation. Click to enlarge.</h4>
 
</div><!-- image8 -->
 
</div><!-- center -->
 
<p>
 
The technique involves defining
 
linear functions connecting consecutive data points to fill in
 
the gaps, using (B3.1).
 
</p>
 
 
<div id="eq17">
 
<p>
 
$$
 
\hat{f}(x)=f(x_i)+\frac{f(x_{i+1})-f(x_i)}{x_{i+1}-x_i}
 
\cdot (x-x_i) \; \text{for} \; x_i < x < x_{i+1} \;\;\; \text{(B3.1)}
 
$$
 
</p>
 
</div>
 
<p>
 
Bilinear interpolation is the generalization
 
of this interpolation technique to two dimensions,
 
illustrated using Figure 9.
 
</p>
 
 
  <div class="center">
 
<div id="image9">
 
<a class="example-image-link"
 
href="https://static.igem.org/mediawiki/2015/c/c1/KU_Leuven_bilinear_interpolation.png"
 
data-lightbox="example-set"
 
data-title="Illustration of bilinear interpolation">
 
<img class="example-image"
 
src="https://static.igem.org/mediawiki/2015/c/c1/KU_Leuven_bilinear_interpolation.png"
 
alt="Illustration of bilinear interpolation"
 
width="60%"
 
></a>
 
<h4>
 
<div id=figure9>Figure 9</div>
 
Illustration of bilinear interpolation. Click to enlarge.</h4>
 
</div><!-- image9 -->
 
</div><!-- center -->
 
<p>
 
 
Let’s assume we want to interpolate
 
a two-dimensional function within a rectangular domain where
 
the function values are only known on the four adjacent grid
 
points $Q$. First we linearly interpolate in the x-direction on
 
the lines connecting $Q_{i,j}$ and $Q_{i+1,j}$,
 
$Q_{i,j+1}$ and $Q_{i+1,j+1}$,
 
which yields (B3.2).
 
</p>
 
 
<div id="eq18">
 
<p>
 
$$
 
\hat{f}(x,y_j)=\frac{x_{i+1}-x}{x_{i+1}-x_i} \cdot f(Q_{i,j})
 
+ \frac{x-x_i}{x_{i+1}-x_i} \cdot f(Q_{i+1,j})
 
\;\;\; \text{(B3.2a)}
 
$$
 
$$
 
\hat{f}(x,y_{j+1})=\frac{x_{i+1}-x}{x_{i+1}-x_i} \cdot f(Q_{i,j+1})
 
+ \frac{x-x_i}{x_{i+1}-x_i} \cdot f(Q_{i+1,j+1})
 
\;\;\; \text{(B3.2b)}
 
$$
 
</p>
 
</div>
 
<p>
 
We then obtain interpolated
 
values for all positions on the bottom and top edge. The key
 
idea then is to take these values and interpolate in the y-direction
 
for any given x-position, yielding (B3.3).
 
</p>
 
 
<div id="eq19">
 
<p>
 
$$
 
\hat{f}(x,y)=\frac{y_{j+1}-y}{y_{j+1}-y_j} \cdot \hat{f}(x,y_j)
 
+ \frac{y-y_j}{y_{j+1}-y_j} \cdot \hat{f}(x,y_{j+1}) \\
 
 
=\frac{y_{j+1}-y}{y_{j+1}-y_j} \cdot \Bigg[
 
\frac{x_{i+1}-x}{x_{i+1}-x_i} \cdot f(Q_{i,j})
 
+ \frac{x-x_i}{x_{i+1}-x_i} \cdot f(Q_{i+1,j})
 
\Bigg] \\
 
 
+\frac{y-y_j}{y_{j+1}-y_j} \cdot \Bigg[
 
\frac{x_{i+1}-x}{x_{i+1}-x_i} \cdot f(Q_{i,j+1})
 
+ \frac{x-x_i}{x_{i+1}-x_i} \cdot f(Q_{i+1,j+1})
 
\Bigg]
 
\;\;\; \text{(B3.3)}
 
$$
 
</p>
 
</div>
 
<p>
 
This then gives
 
interpolated values for any point within the domain. Note that
 
the resulting function is not linear, but quadratic since it contains the term
 
$x \cdot y$.
 
</p>
 
</div><!-- widebox -->
 
</div><!-- toggleseven -->
 
</div><!-- togglebar -->
 
 
</div><!-- center -->
 
 
<div class="part">
 
<h2>Synthesis</h2>
 
<p>
 
The diffusion, production and degradation of AHL and leucine are
 
described by the PDEs (4). At the same time, the random movement,
 
chemotaxis and intercellular
 
interactions of cell types A and B are captured by
 
the stochastic equation of motion (11). Translating the spatial distribution of
 
cells type A into a density field to feed into the PDE module can be
 
accomplished by using kernel density estimation (16). Finally, the agents
 
can request concentrations and gradients at arbitrary positions from the
 
PDE module using
 
bilinear interpolation (B3.3). Taken together, these equations describe the
 
individual modules, as well as the two-way information transfer in between them,
 
and thus they fully define our hybrid model. A graphical summary is shown in
 
Figure 10.
 
 
</p>
 
<div class="whiterow"></div>
 
 
  <div class="center">
 
<div id="image10">
 
<a class="example-image-link"
 
href="https://static.igem.org/mediawiki/2015/b/b3/KU_Leuven_hybrid_model_summary.jpg"
 
data-lightbox="example-set"
 
data-title="Graphical summary of hybrid model">
 
<img class="example-image"
 
src="https://static.igem.org/mediawiki/2015/b/b3/KU_Leuven_hybrid_model_summary.jpg"
 
alt="Graphical summary of hybrid model"
 
width="80%">
 
</a>
 
<h4>
 
<div id=figure10>Figure 10</div>
 
Graphical summary of hybrid model. Click to enlarge.</h4>
 
</div><!-- image10 -->
 
</div><!-- center -->
 
</div><!-- part -->
 
 
<div class="whiterow"></div>
 
</div><!-- summarytext1 -->
 
 
 
<!--</div>-->
 
<!--</div>-->
 
 
<!-- new section -->
 
<div class="summaryheader">
 
<div class="summaryimg">
 
<img src="https://static.igem.org/mediawiki/2015/5/5c/KU_Leuven_Banner_Groen2.jpg"
 
width="100%">
 
<div class="head">
 
<h2>Implementation</h2>
 
</div><!-- head -->
 
</div><!-- summaryimg -->
 
</div><!-- summaryheader -->
 
 
<div class="summarytext1">
 
<div class="part">
 
<h2>Agent-Based Module</h2>
 
<p>
 
Introductory text about the Agent-Based Module
 
</p>
 
</div><!-- part -->
 
 
<div class="center"><!-- start of toggles -->
 
<div class="togglebar">
 
<div class="toggleeight">
 
<h2>Nearest-Neighbor Algorithm</h2>
 
</div>
 
<div id="toggleeight">
 
               
 
<p>
 
The cell-cell interactions have already been fully described in the paragraphs above. However, solving the equation of motion of an agent in its current form requires the computation of the force due to every other agent. If we take N to be the number of agents, that means $N*(N-1)$ amount of force computations are needed in total and therefore the computation time grows as $O(N^2)$. To put this in perspective, if we simulate a thousand agents, the amount of interactions adds up to around one million. This puts a heavy computational load on the model which can be reduced greatly by using a smarter algorithm. The crucial observation here is that the force goes to zero when the distance between two cells is larger than $2\cdot r_{cutoff}$. An agent therefore only interacts with its nearest neighbors, meaning that the naive implementation wastes a lot of time computing interactions which have no effect anyways. Moreover, as the simulation progresses, an agent’s neighbors do not vary greatly from one timestep to another. A more efficient approach then consists of periodically searching for the nearest neighbors within a fixed radius, storing them in a list for every agent and updating their coordinates for all intermediate timesteps.
 
 
Since the agents are programmed to repel each other when they approach each other too closely, they will evolve to a rather uniform distribution. The most appropriate fixed-radius nearest neighbor search algorithm in that case is the cell technique. In this algorithm, the agents are structured by placing a rectangular grid of cells over the domain and assigning every agent to a cell (fig9). Determining which cell an agent belongs to is easily done by rounding down the x and y-coordinates to the nearest multiple of the grid spacing. If the grid spacing is chosen so that interactions do not reach further than adjacent cells, every agent only needs to look for neighbors in 9 cells (its own cells plus 8 adjacent cells) instead of the entire domain.
 
 
 
 
</p>
 
 
</div><!-- toggleeight -->
 
</div><!-- togglebar -->
 
 
<div class="togglebar">
 
<div class="togglenine">
 
<h2>Choice of Timestep</h2>
 
</div>
 
<div id="togglenine">
 
<p>
 
In contrast to PDE models, there is no danger of instability
 
due to a large timestep in agent-based models. However, since
 
the time step determines how far a cell can “jump” away from
 
its previous position, we need to ensure that it is not so
 
large as to negate the effect of intercellular interactions.
 
To understand this, we first discretize the Brownian motion term
 
in the agent’s equation of motion (17).
 
</p>
 
<div id="eq23">
 
<p>
 
$$
 
\sqrt{2 \cdot \mu}\cdot d\vec{W} \rightarrow
 
\sqrt{2 \cdot \mu \cdot \Delta t}\cdot \vec{N}(0,1)
 
\;\;\; \text{(17)}
 
$$
 
</p>
 
</div>
 
<p>
 
The details of this discretization scheme goes beyond the scope of this
 
text and will not be given here.
 
It can be proven that
 
this term follows a distribution as (18).
 
</p>
 
<div id="eq24">
 
<p>
 
$$
 
\sqrt{2 \cdot \mu \cdot \Delta t}\cdot \vec{N}(0,1) \sim
 
\vec{N}(0,2 \cdot \mu \cdot \Delta t)
 
\;\;\; \text{(18)}
 
$$
 
</p>
 
</div>
 
<p>
 
Now assume that there
 
are two A type cells who are exactly $2 \cdot r_0$ apart.
 
If the timestep
 
is too large, they might “jump” out of each other’s influence
 
at the next iteration. This situation is highly unrealistic
 
because the cells actually should feel each other’s pull while
 
they randomly drift away from each other. If the timestep is
 
too large, there’s a high probability that the cells detach
 
without experiencing the potential at all and therefore the
 
probability of detachment is artificially inflated.
 
 
Of course, we can never completely eliminate the possibility
 
of a sudden detachment because the random is sampled by a
 
normal distribution, but we can define a desired probability
 
of sudden detachment $P^-$. Then, for a given $r_0$,
 
$r_{cutoff}$ and $\mu_A$,
 
we can determine how small the timestep has to be for the
 
sudden detachment to happen only with probability $P^-$. First,
 
we assume two A type cells with coordinates $(x_1,y_1)$ and
 
$(x_2,y_2)$ that are at a distance $2 \cdot r_0 = R_0$
 
from each other.
 
We then define the coordinates at the next timestep as
 
stochastic variables $X_1$, $Y_1$, $X_2$ and $Y_2$, which are distributed
 
as (19).
 
</p>
 
<div id="eq25">
 
<p>
 
$$
 
X_1 \sim x_1 + \sqrt{2 \cdot \mu \cdot \Delta t}\cdot {N}(0,1) \sim
 
{N}(x_1,2 \cdot \mu \cdot \Delta t)
 
\;\;\; \text{(19a)}
 
$$
 
$$
 
Y_1 \sim y_1 + \sqrt{2 \cdot \mu \cdot \Delta t}\cdot {N}(0,1) \sim
 
{N}(y_1,2 \cdot \mu \cdot \Delta t)
 
\;\;\; \text{(19b)}
 
$$
 
$$
 
X_2 \sim x_2 + \sqrt{2 \cdot \mu \cdot \Delta t}\cdot {N}(0,1) \sim
 
{N}(x_2,2 \cdot \mu \cdot \Delta t)
 
\;\;\; \text{(19c)}
 
$$
 
$$
 
Y_2 \sim y_2 + \sqrt{2 \cdot \mu \cdot \Delta t}\cdot {N}(0,1) \sim
 
{N}(y_2,2 \cdot \mu \cdot \Delta t)
 
\;\;\; \text{(19d)}
 
$$
 
</p>
 
</div>
 
<p>
 
We are not interested in the future coordinates
 
as such, but rather in the future intercellular distance,
 
which is composed of these future coordinates and therefore
 
also is a stochastic variable. For ease of calculation, we
 
work with the square of this variable, which we name as R²
 
and define as (20).
 
</p>
 
<div id="eq26">
 
<p>
 
$$
 
R^2=(X_1-X_2)^2 + (Y_1-Y_2)^2
 
\;\;\; \text{(20)}
 
$$
 
</p>
 
</div>
 
<p>
 
It is easy to see that $X_1-X_2$ and
 
$Y_1-Y_2$ are both noncentral (mean is not equal to zero)
 
normally distributed variables as well. Therefore $R^2$ obeys
 
a noncentral chi-squared distribution with 2 degrees of freedom.
 
Before proceeding further, we have to make our $R^2$ conform
 
to a standard noncentral chi-squared distribution. First
 
we normalize $X_1-X_2$ and $Y_1-Y_2$. We observe that $X_1-X_2$
 
and
 
$Y_1-Y_2$ are distributed as (21).
 
</p>
 
<div id="eq27">
 
<p>
 
$$
 
\begin{split}
 
X_1-X_2 & \sim N(x_1-x_2, 2 \cdot \mu \cdot \Delta t
 
+ 2 \cdot \mu \cdot \Delta t) \\
 
& \sim N(x_1-x_2, 4 \cdot \mu \cdot \Delta t)
 
\;\;\; \text{(21a)}
 
\end{split}
 
$$
 
$$
 
\begin{split}
 
Y_1-Y_2 & \sim N(y_1-y_2, 2 \cdot \mu \cdot \Delta t
 
+ 2 \cdot \mu \cdot \Delta t) \\
 
& \sim N(y_1-y_2, 4 \cdot \mu \cdot \Delta t)
 
\;\;\; \text{(21b)}
 
\end{split}
 
$$
 
</p>
 
</div>
 
<p>
 
Secondly, because the
 
means of the $X_1-X_2$ and $Y_1-Y_2$ are not equal to zero, we
 
need to account for this by calculating the noncentrality
 
parameter $\lambda$. This is equal to the sum of squares of
 
the means of the normalized $X_1-X_2$ and $Y_1-Y_2$ variables,
 
which of course equals (22), corresponding to a scaled
 
version of the original intercellular distance.
 
</p>
 
<div id="eq28">
 
<p>
 
$$
 
\lambda = (\frac{x_1-x_2}{2 \cdot \sqrt{\mu \cdot \Delta t})^2
 
+ (\frac{y_1-y_2}{2 \cdot \sqrt{\mu \cdot \Delta t})^2 \\
 
\frac{R^2_0}{4 \cdot \mu \Delta t}
 
\;\;\; \text{(22)}
 
$$
 
</p>
 
</div>
 
<p>
 
We then
 
define the standardized variable R’² as (eq. 21). To recap,
 
the meaning of R’² is the scaled squared intercellular distance
 
after updating once with timestep \triangle t starting from
 
an initial distance of R0. Then we define F as the cumulative
 
distribution function (eq. 22). To then finally calculate
 
appropriate timestep for a given r0, rcutoff and µA we take
 
the inverse cumulative distribution function and equate it
 
to the scaled squared intercellular cutoff distance (eq. 23).
 
Note that both the left hand side and right hand side
 
depend on the \triangle t in a nonlinear way, therefore
 
a nonlinear solver is required.
 
</p>
 
</div><!-- togglenine -->
 
</div><!-- togglebar -->
 
</div><!-- center, end of toggles -->
 
 
<div class="part">
 
<h2>Partial Differential Equations Module</h2>
 
<p>
 
Introductory text about the Partial Differential Equations Module
 
</p>
 
</div><!-- part -->
 
 
<div class="center">
 
<div class="togglebar">
 
<div class="toggleten">
 
<h2>Computational Scheme</h2>
 
</div><!-- toggleten -->
 
<div id="toggleten">
 
                <p>
 
                As mentioned earlier the concentrations of AHL and Leucine are modeled using partial differential equations.
 
                In the colony level model these equations are solved explicitly. Explicit schemes do not require a lot of
 
                work per time step, but unfortunately are not unconditionally stable. In two dimensions the grid ratios
 
                $dt/dx^2$ and $dt/dy^2$ can not exceed $dt/dx^2 + dt/dy^2 \leq \frac{1}{2}$ for the solver to be stable.
 
                When computing the solution of the hybrid model this requrement forces us to spend a lot of CPU time solving
 
                partial differential equations that could be better spent simulation the agents. Therefore an implicit ADI
 
                Alternating direction implicit scheme has been implemented.  ADI-schemes are unconditionally stable, which
 
                allows it to take large time steps with the PDE solver. We used the following scheme:
 
                          $$ (1 - \frac{1}{2} \mu_x \delta_x^2) U^{n+\frac{1}{2}} + \frac{1}{4}kU^{n+\frac{1}{2}}
 
                          = (1 + \frac{1}{2} \mu_y \delta_y^2) U^n - \frac{1}{4}kU^n + \frac{\alpha}{2} \rho_A  $$
 
                          $$ (1 - \frac{1}{2} \mu_y \delta_y^2) U^{n+1} + \frac{1}{4}kU^{n+1} =
 
                          (1 + \frac{1}{2} \mu_x \delta_x^2)U^{n+\frac{1}{2}} - \frac{1}{4}kU^{n+\frac{1}{2}} + \frac{\alpha}{2} \rho_A $$
 
                In the equations above $\mu$ denotes grid ratios and $\delta^2$ central differences. The production and
 
                degradation terms have been incorporated at every time level with a factor of $\frac{1}{4}$. 
 
                The image below shows the computational molecule of the ADI scheme we chose to implement:
 
                <br/>
 
                </p>
 
 
<div class="center">
 
<div id="image2">
 
                        <a class="example-image-link" href="https://static.igem.org/mediawiki/2015/f/f7/KU_Leuven_ADI_Molecule.png" data-lightbox="example-set" data-title="Epanechnikov kernel with h=1"><img class="example-image" src="https://static.igem.org/mediawiki/2015/f/f7/KU_Leuven_ADI_Molecule.png" alt="Epanechnikov kernel with h=1" width="45%" height="45%"></a>
 
                        <h4><div id=figure1>Figure 4</div> ADI-Molecule. Click to enlarge.</h4>
 
</div><!-- image2 -->
 
</div><!-- center -->
 
 
</div><!-- toggleten -->
 
</div><!-- togglebar -->
 
 
<div class="togglebar">
 
<div class="toggleeleven">
 
<h2>Boundary Conditions</h2>
 
</div><!-- toggleeleven -->
 
<div id="toggleeleven">
 
<p>
 
Text about boundary conditions.
 
</p>
 
</div><!-- toggleeleven -->
 
</div><!-- togglebar -->
 
</div><!-- center -->
 
 
 
<div class="part">
 
<h2>Matching</h2>
 
<p>
 
Introductory text about matching
 
</p>
 
</div><!-- part -->
 
 
 
<div class="center"><!-- start of toggles -->
 
<div class="togglebar">
 
<div class="toggletwelve">
 
            <h2>Decoupling of Timesteps</h2>
 
</div>
 
<div id="toggletwelve">
 
              <p>
 
              In order to benefit from the implicit PDE-solver described above the agent's time steps are chosen smaller
 
              then the time steps of the PDE solver. However type A cells produce molecules continuously as they move
 
              trough space. Therefore we record their positions and average over the kernel functions centered at
 
              past positions since the last PDE evaluation. That way we avoid blurring the results of the PDE solver
 
              too much, when the time step of the agents is reduced.
 
              Last but not least we want to point out the relationship between kernel bandwidth and the grid on which the
 
              PDE is solved. The larger the kernel bandwidth $h$ is chosen the coarser the PDE grid can be without loosing
 
              cells between grid points. If the the diameter of the kernel functions is smaller then the distance between
 
              PDE grid points it can happen, that the kernel does not overlap with any PDE grid points. In this case a cells
 
              contribution to the molecule concentrations at the next time step could be lost. When the PDE grid is
 
              widened to save computing time it is thus necessary to increase the Kernel bandwidth as well. However this
 
              increase can lead to a situation where the Kernel function covers significantly more space then the actual
 
              diameter of a bacterium. In these cases the Kernel function can be interpreted as a probability function for
 
              a cells position. However we avoided too large bandwiths and PDE grids by running our 2-D simulations at the
 
              Flemish Supercomputer Center (VSC).
 
              </p>
 
 
 
                  <div class="center">
 
                        <div id="image2">
 
                        <a class="example-image-link" href="https://www.vscentrum.be/img/svg/logo.svg" data-lightbox="example-set" data-title="Logo of the Flemish Supercomputing Center (VSC)"><img class="example-image" src="https://www.vscentrum.be/img/svg/logo.svg" alt="Logo Flemish Supercomputer Center" width="30%" height="35%"></a>
 
                        <h4><div id=figure3>Figure 6</div> Logo of the Flemish Supercomputer Center. Click to enlarge.</h4>
 
                        </div>
 
</div><!-- center -->
 
</div><!-- toggletwelve -->
 
</div><!-- togglebar -->
 
 
<div class="togglebar">
 
<div class="toggle13">
 
<h2>Bandwidth Tuning</h2>
 
</div>
 
<div id="toggletwelve">
 
<p>
 
Text about bandwidth tuning.
 
</p>
 
</div><!-- toggle13 -->
 
</div><!-- togglebar -->
 
</div><!-- center, end of toggles -->
 
 
<div class="whiterow"></div>
 
</div>
 
<div class="summaryheader">
 
                <div class="summaryimg">
 
                    <img src="https://static.igem.org/mediawiki/2015/5/5c/KU_Leuven_Banner_Groen2.jpg"
 
                        width="100%">
 
                    <div class="head">
 
    <h2>1-D Hybrid Model</h2>
 
</div>
 
</div>
 
</div>
 
<div class="summarytext1">
 
<div class="part">
 
 
            <!-- first Videobox start-->
 
            <video id="video1" preload="auto" width="65%" tabindex="0" controls="" type="video/mp4">
 
              <source type="video/mp4" src="https://static.igem.org/mediawiki/2015/d/df/KU_Leuven_hybrid_1-D.mp4">
 
              Sorry, your browser does not support HTML5 video.
 
            </video>
 
            </br>
 
            <button type="button" onclick="Set1()">Constant speed</button>
 
            <button type="button" onclick="Set2()">Random step</button>
 
            <script>
 
            function Set1() {
 
 
            document.querySelector("#video1 > source").src = "https://static.igem.org/mediawiki/2015/d/df/KU_Leuven_hybrid_1-D.mp4"
 
            document.querySelector("#video1").load();
 
            document.querySelector("#video1").play();
 
            }
 
            function Set2() {
 
 
            document.querySelector("#video1 > source").src = "https://static.igem.org/mediawiki/2015/3/3d/KU_Leuven_1-
 
            D_random_step.mp4"
 
            document.querySelector("#video1").load();
 
            document.querySelector("#video1").play();
 
            }     
 
            </script>
 
            <!-- video end -->
 
            <p> </br>The video box above shows one dimensional simulation results for the hybrid model. A constant speed and random step
 
                simulation has been computed. We observe that the bacteria form a traveling wave in both cases, which is essential
 
                for pattern formation. These results are also similar to what we get from the continuous  model, which confirms our
 
                results.
 
            </p>
 
 
</div>
 
</div>
 
<div class="summaryheader">
 
                <div class="summaryimg">
 
                    <img src="https://static.igem.org/mediawiki/2015/5/5c/KU_Leuven_Banner_Groen2.jpg"
 
                        width="100%">
 
                    <div class="head">
 
    <h2>2-D Hybrid Model</h2>
 
</div>
 
</div>
 
</div>
 
<div class="summarytext1">
 
<div class="part">
 
 
            <!-- second Videobox start-->
 
            <video id="video2" preload="auto" width="65%" tabindex="0" controls="" type="video/mp4">
 
            <source type="video/mp4" src="https://static.igem.org/mediawiki/2015/4/4d/KU_Leuven_Hybrid5_2d_clean.mp4">
 
            Sorry, your browser does not support HTML5 vido.
 
            </video>
 
            </br>
 
            <button type="button" onclick="Set3()">Initial condition 1</button>
 
            <button type="button" onclick="Set4()">Initial condition 2</button>
 
            <button type="button" onclick="Set5()">Initial condition 1</button>
 
            <button type="button" onclick="Set6()">Initial condition 2</button>
 
            <button type="button" onclick="Set7()">Random Data</button>
 
            <script>
 
            function Set3() {
 
 
            document.querySelector("#video2 > source").src = "https://static.igem.org/mediawiki/2015/b/b0/KU_Leuven_9VSC.mp4"
 
            document.querySelector("#video2").load();
 
            document.querySelector("#video2").play();
 
            }
 
            function Set4() {
 
 
            document.querySelector("#video2 > source").src = "https://static.igem.org/mediawiki/2015/9/94/KU_Leuven_StartVSC.mp4"
 
            document.querySelector("#video2").load();
 
            document.querySelector("#video2").play();
 
            }
 
            function Set5() {
 
 
            document.querySelector("#video2 > source").src = "https://static.igem.org/mediawiki/2015/b/bf/KU_Leuven_9VSC2.mp4"
 
            document.querySelector("#video2").load();
 
            document.querySelector("#video2").play();
 
            }
 
            function Set6() {
 
 
            document.querySelector("#video2 > source").src = "https://static.igem.org/mediawiki/2015/e/e5/KU_Leuven_StarVSC2.mp4"
 
            document.querySelector("#video2").load();
 
            document.querySelector("#video2").play();
 
            }
 
            function Set7() {
 
 
            document.querySelector("#video2 > source").src = "https://static.igem.org/mediawiki/2015/4/4d/KU_Leuven_Hybrid5_2d_clean.mp4"
 
            document.querySelector("#video2").load();
 
            document.querySelector("#video2").play();
 
            }     
 
            </script>
 
            <!-- video end -->
 
            <p></br> The videos above show simulation videos computed at the Flemish supercomputing center, for three different
 
            initial conditions similar to the ones we used for the colony level model. The first and second condition start
 
            from 9 mixed or 5 colonies of both cell types, arranged in a block or star shape. These first two gradually
 
            separate in a manner similar to what we would we also saw in the colony level model. The result for random
 
            initial data is fundamentally different. As the agent based approach allows for better implementation
 
            of adhesion large cell type A bands form. The AHL and Leucine produced by the type A bacteria causes the
 
            B type cells to move away leading to a pattern which we could not produce using PDEs alone, this
 
            beautifully illustrates the added value of hybrid modeling.</p>
 
      </div>
 
</div>
 
 
<div class="summaryheader">
 
                <div class="summaryimg">
 
                    <img src="https://static.igem.org/mediawiki/2015/5/5c/KU_Leuven_Banner_Groen2.jpg"
 
                        width="100%">
 
                    <div class="head">
 
    <h2>Incorporation of internal model</h2>
 
</div>
 
</div>
 
</div>
 
<div class="summarytext1">
 
<div class="part">
 
      <p>Up until now, we have largely ignored the inner life of the bacteria. This inner life consists of transcriptional  networks and protein kinetics. Instead we assumed that AHL and leucine production is directly proportional to the density of type A cells. This only works in theory, since bacteria will be affected by their surroundings and the way their dynamics react to it.  For example bacteria surrounded by a large concentration of AHL, will have more CheZ and will react more on the presence of Leucine. Also bacteria have different histories and will have different levels of transcription factors and different levels of proteins in their plasma. The proteins are not directly degraded and will still be present in the cytoplasm of the bacteria long after the network has been deactivated. From this, it is clear that 2 bacteria, although surrounded by the same AHL and leucine concentrations, can show different behavior and reaction kinetics. </p>
 
<br>
 
 
      <p>This results in a heterogeneity of the bacterial population that has not yet been accounted for. To make up for this anomaly, we decided to add an internal model to every agent. This way we will get more realistic simulations. Every agent will get their own levels of CheZ, LuxR, LuxI and so on and will have individual reactions on their surroundings. We hope that this way we can get closer to the behavior of real bacteria. </p>
 
</div>
 
</div>
 
 
<!--References-->
 
<div class="summaryheader">
 
    <div class="summaryimg">
 
  <img src="https://static.igem.org/mediawiki/2015/5/5c/KU_Leuven_Banner_Groen2.jpg" width="100%">
 
  <div class="head">
 
      <h2> References </h2>
 
    </div>
 
  </div>
 
  </div>
 
<div class="reference">
 
<div class="part">
 
       
 
      <table>
 
      <tr valign="top">
 
        <td align="right" class="bibtexnumber">
 
        [<a name="Franz2013">1</a>]
 
        </td>
 
        <td class="bibtexitem">
 
        Benjamin Franz and Radek Erban.
 
        Hybrid modelling of individual movement and collective behaviour.
 
          <em>Lecture Notes in Mathematics</em>, 2071:129-157, 2013.
 
          [&nbsp;<a href="http://link.springer.com/chapter/10.1007%2F978-3-642-35497-7_5">.pdf</a>&nbsp;]
 
 
          </td>
 
          </tr>
 
 
 
          <tr valign="top">
 
          <td align="right" class="bibtexnumber">
 
          [<a name="Guo2008">2</a>]
 
          </td>
 
          <td class="bibtexitem">
 
          Zaiyi Guo, Peter M&nbsp;A Sloot, and Joc&nbsp;Cing Tay.
 
          A hybrid agent-based approach for modeling microbiological systems.
 
          <em>Journal of Theoretical Biology</em>, 255(2):163-175, 2008.
 
          [&nbsp;<a href="http://dx.doi.org/10.1016/j.jtbi.2008.08.008">DOI</a>&nbsp;]
 
 
          </td>
 
          </tr>
 
 
 
          <tr valign="top">
 
          <td align="right" class="bibtexnumber">
 
          [<a name="Keller1971">3</a>]
 
          </td>
 
          <td class="bibtexitem">
 
          E&nbsp;F Keller and L&nbsp;A Segel.
 
          Traveling bands of chemotactic bacteria: a theoretical analysis.
 
          <em>Journal of theoretical biology</em>, 30(2):235-248, 1971.
 
          [&nbsp;<a href="http://dx.doi.org/10.1016/0022-5193(71)90051-8">DOI</a>&nbsp;]
 
 
          </td>
 
          </tr>
 
 
 
          <tr valign="top">
 
          <td align="right" class="bibtexnumber">
 
          [<a name="Purcell1977">4</a>]
 
          </td>
 
          <td class="bibtexitem">
 
          E.&nbsp;M. Purcell.
 
          Life at low Reynolds number, 1977.
 
          [&nbsp;<a href="http://dx.doi.org/10.1119/1.10903">DOI</a>&nbsp;]
 
 
          </td>
 
          </tr>
 
 
 
          <tr valign="top">
 
          <td align="right" class="bibtexnumber">
 
          [<a name="Stevens2000">5</a>]
 
          </td>
 
          <td class="bibtexitem">
 
          Angela Stevens.
 
          The Derivation of Chemotaxis Equations as Limit Dynamics of
 
            Moderately Interacting Stochastic Many-Particle Systems, 2000.
 
          [&nbsp;<a href="http://dx.doi.org/10.1137/S0036139998342065">DOI</a>&nbsp;]
 
 
          </td>
 
          </tr>
 
 
          <tr valign="top">
 
          <td align="right" class="bibtexnumber">
 
          [<a name="Schaefer2000">6</a>]
 
          </td>
 
          <td class="bibtexitem">
 
    A.&nbsp;L. Schaefer, B.&nbsp;L. Hanzelka, M.&nbsp;R. Parsek, and E.&nbsp;P. Greenberg.
 
    Detection, purification, and structural elucidation of the
 
    acylhomoserine lactone inducer of Vibrio fischeri luminescence and other
 
    related molecules.
 
    <em>Bioluminescence and Chemiluminescence, Pt C</em>, 305:288-301,
 
    2000.
 
          </td>
 
          </tr>
 
      </table>
 
    </div>
 
</div>
 
 
 
<div class="summarytext1">
 
<div class="part">
 
</div>
 
</div>
 
 
<div class="whiterow"></div>
 
            <!--Foot don't touch!-->
 
 
            <div class="subsections">
 
                    <div class="subsectionwrapper">
 
                        <div class="subimgrow">
 
                         
 
                            <div class="subimg">
 
                                <a href="https://2015.igem.org/Team:KU_Leuven/Modeling/Top">
 
                                    <img
 
                                        src="https://static.igem.org/mediawiki/2015/6/6a/KU_Leuven_Wiki_Button_-_Colony_level2.png"
 
                                        width="100%"></a>
 
                            </div>
 
 
                            <div class="whitespace"></div>
 
 
                            <div class="subimg">
 
                                <a href="https://2015.igem.org/Team:KU_Leuven/Modeling/Internal">
 
                                    <img
 
                                        src="https://static.igem.org/mediawiki/2015/4/47/KU_Leuven_Wiki_Button_-_Internal_model2.png"
 
                                        width="100%"></a>
 
                            </div>
 
 
                            <div class="whitespace"></div>
 
 
                            <div class="subimg">
 
                                <a href="https://2015.igem.org/Team:KU_Leuven/Modeling/Toulouse">
 
                                    <img src="https://static.igem.org/mediawiki/2015/b/b1/KU_Leuven_Wiki_Button_Flux2.png"
 
                                        width="100%"></a>
 
                            </div>
 
 
                            <div class="whitespace"></div>
 
 
                            <div class="subimg">
 
        <a href="https://2015.igem.org/Team:KU_Leuven/Modeling" >
 
    <img src="https://static.igem.org/mediawiki/2015/c/cb/KUL_Wiki_Button_-_Back.png" width="100%" >
 
                                </a>
 
    </div>
 
 
                           
 
                        </div>
 
 
                        <div class="subtextrow">
 
                         
 
                            <div class="subtext">
 
                                <a href="https://2015.igem.org/Team:KU_Leuven/Modeling/Top">
 
                                    <h2>Colony</h2>
 
                                    <p>
 
                                        Our colony layer model relies on a Keller-Segel type system of differential
 
                                        equations. These equations are simulated using finite differences.
 
                                     
 
                                    </p>
 
                                </a>
 
                            </div>
 
 
                            <div class="whitespace"></div>
 
 
                            <div class="subtext">
 
                                <a href="https://2015.igem.org/Team:KU_Leuven/Modeling/Internal">
 
                                    <h2>Internal</h2>
 
                                    <p>Our internal model aims to simulate the internal dynamics of every cell with a system of ordinary differential equations.</p>
 
                                </a>
 
                            </div>
 
 
                            <div class="whitespace"></div>
 
 
                            <div class="subtext">
 
                                <a href="https://2015.igem.org/Team:KU_Leuven/Modeling/Toulouse">
 
                                    <h2>FBA (Toulouse)</h2>
 
                                    <p>
 
                                        As a part of our modeling cooperation we exchanged models with the Toulouse
 
                                        team. This is the flux balance analysis they performed for us.<br/>
 
                                    </p>
 
                                </a>
 
                            </div>
 
 
                            <div class="whitespace"></div>
 
 
                            <div class="subtext">
 
                                <a href = "https://2015.igem.org/Team:KU_Leuven/Modeling">
 
            <h2>Back</h2>
 
    <p>
 
        Go back to the Modeling page.
 
    </p>
 
</a>
 
                            </div>
 
                         
 
                        </div>
 
 
                        <div class="subimgreadmore">
 
                         
 
                            <div class="subimgrm">
 
                                <a href="https://2015.igem.org/Team:KU_Leuven/Modeling/Top">
 
                                    <div id="more">
 
                                        <img alt="Read more" height="40%"
 
                                            src="https://static.igem.org/mediawiki/2015/7/73/KUL_Wiki_Button_-_Read_more.png"
 
                                            width="85%">
 
                                    </div>
 
                                </a>
 
                            </div>
 
 
                            <div class="whitespace"></div>
 
 
                            <div class="subimgrm">
 
                                <a href="https://2015.igem.org/Team:KU_Leuven/Modeling/Internal">
 
                                    <div id="more">
 
                                        <img alt="Read more" height="40%"
 
                                            src="https://static.igem.org/mediawiki/2015/7/73/KUL_Wiki_Button_-_Read_more.png"
 
                                            width="85%">
 
                                    </div>
 
                                </a>
 
                            </div>
 
                            <div class="whitespace"></div>
 
 
                            <div class="subimgrm">
 
                                <a href="https://2015.igem.org/Team:KU_Leuven/Modeling/Toulouse">
 
                                    <div id="more">
 
                                        <img alt="Read more" height="40%"
 
                                            src="https://static.igem.org/mediawiki/2015/7/73/KUL_Wiki_Button_-_Read_more.png"
 
                                            width="85%">
 
                                    </div>
 
                                </a>
 
                               
 
                            </div>
 
 
                            <div class="whitespace"></div>
 
 
                            <div class="subimgrm">
 
                                <a href="https://2015.igem.org/Team:KU_Leuven/Modeling">
 
                                    <div id="back">
 
<img src="https://static.igem.org/mediawiki/2015/7/73/KUL_Wiki_Button_-_Read_more.png" 
 
                                        height="40%" width="85%" alt="Read more">
 
    </div>
 
                                </a>
 
                               
 
                            </div>
 
 
                        </div>
 
                    </div>
 
 
<div class="subimgrowm">
 
                        <div class="whiterow"></div>
 
                    </div>
 
 
<div class="subimgrowm">
 
                        <div class="subimgm">
 
                            <a href="https://2015.igem.org/Team:KU_Leuven/Modeling/Top">
 
                                <b>Colony</b>
 
                                <img
 
                                    src="https://static.igem.org/mediawiki/2015/6/6a/KU_Leuven_Wiki_Button_-_Colony_level2.png"
 
                                    width="100%"></a>
 
                        </div>
 
 
                        <div class="whitespace"></div>
 
 
                        <div class="subtextm">
 
                            <a href="https://2015.igem.org/Team:KU_Leuven/Modeling/Top">
 
                                <p>
 
                                    Our colony layer model relies on a Keller-Segel type system of differential
 
                                    equations. These equations are simulated using finite differences.
 
                                    <br/>
 
                                </p>
 
                            </a>
 
                        </div>
 
                    </div>
 
 
                    <div class="subimgrowm">
 
                        <div class="whiterow"></div>
 
                    </div>
 
 
 
                    <div class="subimgrowm">
 
                        <div class="subimgm">
 
                            <a href="https://2015.igem.org/Team:KU_Leuven/Modeling/Internal">
 
                            <b>Internal</b>
 
                            <img
 
                                src="https://static.igem.org/mediawiki/2015/4/47/KU_Leuven_Wiki_Button_-_Internal_model2.png"
 
                                width="100%"></a>
 
                        </div>
 
 
                        <div class="whitespace"></div>
 
 
                        <div class="subtextm">
 
                            <a href="https://2015.igem.org/Team:KU_Leuven/Modeling/Internal">
 
                                <p>Our internal model aims to simulate the internal dynamics of every cell with a system of ordinary differential equations.
 
                                </p>
 
                            </a>
 
                        </div>
 
                    </div>
 
 
                    <div class="subimgrowm">
 
                        <div class="whiterow"></div>
 
                    </div>
 
 
                  <div class="subimgrowm">
 
 
                        <div class="subimgm">
 
                            <a href="https://2015.igem.org/Team:KU_Leuven/Modeling/Toulouse">
 
                            <b>FBA (Toulouse)</b>
 
                            <img src="https://static.igem.org/mediawiki/2015/b/b1/KU_Leuven_Wiki_Button_Flux2.png"
 
                                width="100%">
 
                            </a>
 
                        </div>
 
 
                        <div class="whitespace"></div>
 
 
                        <div class="subtextm">
 
                              <a href="https://2015.igem.org/Team:KU_Leuven/Modeling/Toulouse">
 
                                <p>
 
                                    As a part of our modeling cooperation we exchanged models with the Toulouse
 
                                        team. This is the flux balance analysis they performed for us.<br/>
 
                                </p>
 
                              </a>                           
 
                        </div>
 
                       
 
                    </div>
 
 
<div class="subimgrowm">
 
<div class="whiterow">
 
</div>
 
</div>
 
 
<div class="subimgrowm">
 
<div class="subimgm">
 
<a href="https://2015.igem.org/Team:KU_Leuven/Modeling">
 
<b>Back</b>
 
<img src="https://static.igem.org/mediawiki/2015/c/cb/KUL_Wiki_Button_-_Back.png" width="100%" >
 
</a>
 
</div>
 
 
<div class="whitespace"></div>
 
 
<div class="subtextm">
 
<a href="https://2015.igem.org/Team:KU_Leuven/Modeling">
 
<p>Go back to the Modeling page.<br/></p>
 
</a>
 
</div>
 
</div>
 
</div>
 
</div>
 
   
 
 
    <div id="summarywrapper">
 
        <div class="summary">
 
            <h3>
 
                Contact
 
            </h3>
 
            <p style="font-size:1.3em; text-align: center">
 
                Address: Celestijnenlaan 200G room 00.08 - 3001 Heverlee<br>
 
                Telephone: +32(0)16 32 73 19<br>
 
                Email: igem@chem.kuleuven.be<br>
 
            </p>
 
        </div>
 
    </div>
 
 
<div class="whiterow"></div>
 
 
    <div id="footer">
 
        <div class="rowfooter">
 
            <div class="emptyfooter"></div>
 
            <div id="bioscenter">
 
                <a href="http://www.kuleuven.be/bioscenter/"><img alt="bioSCENTer"
 
                    src="https://static.igem.org/mediawiki/2015/3/3c/KU_Leuven_Logo_BioSCENTer.png"
 
                    width="95%"></a>
 
            </div>
 
            <div class="emptyfooter"></div>
 
        </div>
 
 
        <div class="rowfooter">
 
            <div id="lrd">
 
                <a href="http://lrd.kuleuven.be/"><img alt="KU Leuven Research & Development"
 
                    src="https://static.igem.org/mediawiki/2015/b/b4/KU_Leuven_HR_LRD_CMYK_LOGO.jpg"
 
                    width="95%"></a>
 
            </div>
 
            <div id="youreca">
 
                <a href="http://www.kuleuven.be/onderzoek/youreca/"><img alt="YOURECA"
 
                    src="https://static.igem.org/mediawiki/2015/4/4f/KU_Leuven_YOURECA.png" width="95%"></a>
 
            </div>
 
            <div id="solvay">
 
                <a href="http://www.solvay.com/"><img alt="SOLVAY"
 
                    src="https://static.igem.org/mediawiki/2015/8/89/KU_Leuven_Solvay_Transparant.png"
 
                    width="95%"></a>
 
            </div>
 
            <div id="medicine">
 
                <a href="http://gbiomed.kuleuven.be/english"><img alt="Faculty of Medicine"
 
                    src="https://static.igem.org/mediawiki/2015/8/83/KU_Leuven_Biomedische_Wetenschappen_logo_.png"
 
                    width="95%"></a>
 
            </div>
 
        </div>
 
 
        <div class="rowfooter">
 
            <div id="imec">
 
                <a href="http://www2.imec.be/be_en/home.html"><img alt="imec"
 
                    src="https://static.igem.org/mediawiki/2015/7/70/KU_Leuven_Imec_logo_transparant.png"
 
                    width="95%"></a>
 
            </div>
 
            <div id="genzyme">
 
                <a href="http://www.genzyme.be/"><img alt="Genzyme"
 
                    src="https://static.igem.org/mediawiki/2015/b/b7/KU_Leuven_Genzyme_logo_transparant.png"
 
                    width="95%"></a>
 
            </div>
 
            <div id="eppendorf">
 
                <a href="https://www.eppendorf.com/BE-en/"><img alt="Eppendorf"
 
                    src="https://static.igem.org/mediawiki/2015/9/96/KU_Leuven_Logo_Eppendorf_transparant.png"
 
                    width="95%"></a>
 
            </div>
 
<div id="saillart">
 
      <a href="http://www.glasatelier-saillart.be/English/english.html"><img src="https://static.igem.org/mediawiki/2015/c/ce/KU_Leuven_Sponsor_Saillard.png" alt="Glasatelier Saillart" width="95%"></a>
 
  </div>
 
            <div id="kuleuven">
 
                <a href="http://www.kuleuven.be/english"><img alt="bioSCENTer"
 
                    src="https://static.igem.org/mediawiki/2015/0/08/KULEUVEN_groot.png" width="95%"></a>
 
            </div>
 
            <div class="spacefooter"></div>
 
            <div id="gips">
 
                <a href="http://www.gipsmineral.de/"><img alt="Gips Mineral"
 
                    src="https://static.igem.org/mediawiki/2015/3/3d/KU_Leuven_Gips_Mineral_logo_transparant.png"
 
                    width="95%"></a>
 
            </div>
 
            <div id="kolo">
 
                <a><img alt="Ko-Lo Instruments"
 
                    src="https://static.igem.org/mediawiki/2015/1/15/KUL_Ko-Lo_Instruments_logo_transparant.png"
 
                    width="95%"></a>
 
            </div>
 
        </div>
 
 
        <div class="rowfooter">
 
            <div id="bioke">
 
                <a href="https://www.bioke.com/"><img alt="Bioké"
 
                    src="https://static.igem.org/mediawiki/2015/e/e1/KUL_Biok%C3%A9_logo_transparant.png"
 
                    width="95%"></a>
 
            </div>
 
<div class="logonormal">
 
<div id="regensys">
 
      <a href="http://regenesys.eu/"><img src="https://static.igem.org/mediawiki/2015/e/eb/KU_Leuven_Logo_Regenesys_Transparant.png" alt="Regenesys" width="95%"></a>
 
  </div>
 
<div class="whiterow"></div>
 
  <div id="thermofisher">
 
      <a href="https://www.fishersci.com/us/en/home.html"><img src="https://static.igem.org/mediawiki/2015/a/aa/KUL_Fischer_Scientific_logo_transparant.png" alt="Thermo Fisher Scientific" width="95%"></a>
 
  </div>
 
</div>
 
<div class="logonormal2">
 
<div id="vwr">
 
      <a href="https://be.vwr.com/store/?&_requestid=866148&_DARGS=/store/cms/be.vwr.com/nl_BE/header_20159241139103.jsp.1_AF&_dynSessConf=4047468000326453053&targetURL=/store/%3F%26_requestid%3D866148&lastLanguage=en&/vwr/userprofiling/EditPersonalInfoFormHandler.updateLocale=&_D%3AcurrentLanguage=+&currentLanguage=en&_D%3AlastLanguage=+&_D%3A/vwr/userprofiling/EditPersonalInfoFormHandler.updateLocale=+"><img src="https://static.igem.org/mediawiki/2015/8/8d/KU_Leuven_Logo_VWR_transparant_.png" alt="VWR" width="95%"></a>
 
  </div>
 
<div class = "whiterow"></div>
 
<div id="lgc">
 
<a href="http://www.lgcgroup.com/our-science/genomics-solutions/#.Vfx9V9yLTIU">
 
                <img src="https://static.igem.org/mediawiki/2015/e/e6/KU_Leuven_LOGO_LGC.png" alt="LGC Genomics" width="80%">
 
</a>
 
</div>
 
</div>
 
            <div id="footerimg">
 
                <img
 
                    src="https://static.igem.org/mediawiki/2015/b/b9/KU_Leuven_Zebra_spots_wiki_footer_main.png"
 
                    width="95%">
 
            </div>
 
        <div class="logonormal2">
 
<div id="gimv">
 
      <a href="http://www.gimv.com/en"><img src="https://static.igem.org/mediawiki/2015/a/ac/KU_Leuven_Logo_Gimv_Transparant.png" alt="Gimv" width="95%"></a>
 
  </div>
 
<div class = "whiterow"></div>
 
<div id="sopach">
 
      <a href="http://www.sopachem.com/"><img src="https://static.igem.org/mediawiki/2015/5/55/KU_Leuven_Sopachem.jpeg" alt="Sopachem" width="95%"></a>
 
  </div>
 
</div>
 
            <div id="machery">
 
                <a href="http://www.filterservice.be/"><img alt="Machery Nagel"
 
                    src="https://static.igem.org/mediawiki/2015/4/41/KU_Leuven_Macherey_Nagel_logo_transparant.png"
 
                    width="95%"></a>
 
            </div>
 
<div class="logosmall">
 
    <div id="sigma">
 
      <a href="https://www.sigmaaldrich.com/belgium-nederlands.html"><img src="https://static.igem.org/mediawiki/2015/4/4b/KUL_Sigma-Aldrich_logo_transparant.png" alt="Sigma-Aldrich" width="95%"></a>
 
    </div>
 
    <div class="whiterow"></div>
 
    <div id="egilab">
 
      <a href="http://www.egilabo.be/"><img src="https://static.igem.org/mediawiki/2015/e/e9/KUL_Egilabo_logo_transparant.png" alt="Egilabo" width="95%"></a>
 
    </div>
 
    <div class="whiterow"></div>
 
    <div id="novolab">
 
      <a href="https://www.novolab.be/"><img src="https://static.igem.org/mediawiki/2015/4/4c/KU_Leuven_Novalab.png" alt="Novolab" height="95%"></a>
 
    </div>
 
  </div>
 
        </div>
 
</body>
 
 
<script>
 
$("document").ready(function(){
 
$("#toggleone").hide();
 
$("#toggletwo").hide();
 
$("#togglethree").hide();
 
$("#togglefour").hide();
 
$("#togglefive").hide();
 
$("#togglesix").hide();
 
$("#toggleseven").hide();
 
$("#toggleeight").hide();
 
$("#togglenine").hide();
 
$("#toggleten").hide();
 
$("#toggleeleven").hide();
 
$("#toggletwelve").hide();
 
$("#toggle13").hide();
 
$("#toggle14").hide();
 
 
 
});
 
 
</script>
 
<script>
 
$(".toggleone").click(function () {
 
    $("#toggleone").toggle();
 
});
 
$(".toggletwo").click(function () {
 
    $("#toggletwo").toggle();
 
});
 
$(".togglethree").click(function () {
 
    $("#togglethree").toggle();
 
});
 
$(".togglefour").click(function () {
 
    $("#togglefour").toggle();
 
});
 
$(".togglefive").click(function () {
 
    $("#togglefive").toggle();
 
});
 
$(".togglesix").click(function () {
 
    $("#togglesix").toggle();
 
});
 
$(".toggleseven").click(function () {
 
    $("#toggleseven").toggle();
 
});
 
$(".toggleeight").click(function () {
 
    $("#toggleeight").toggle();
 
});
 
$(".togglenine").click(function () {
 
    $("#togglenine").toggle();
 
});
 
$(".toggleten").click(function () {
 
    $("#toggleten").toggle();
 
});
 
$(".toggleeleven").click(function () {
 
    $("#toggleeleven").toggle();
 
});
 
$(".toggletwelve").click(function () {
 
    $("#toggletwelve").toggle();
 
});
 
$(".toggle13").click(function () {
 
    $("#toggle13").toggle();
 
});
 
$(".toggle14").click(function () {
 
    $("#toggle14").toggle();
 
});
 
</script>
 
 
<!-- sidebar scripts -->
 
 
<script>
 
$(document).ready(function(){
 
 
$(".sidebarright").hide();
 
$(".sidebarleft").hide();
 
 
});
 
 
$( "#eq1" ).mouseenter(function() {
 
$(".sidebarright").show();
 
$("#eq1text").show();
 
$("#eq2text").hide();
 
$("#eq3text").hide();
 
$("#eq4text").hide();
 
$("#eq5text").hide();
 
$("#eq6text").hide();
 
$("#eq7text").hide();
 
$("#eq8text").hide();
 
$("#eq9text").hide();
 
$("#eq10text").hide();
 
$("#eq11text").hide();
 
$("#eq12text").hide();
 
$("#eq13text").hide();
 
$("#eq14text").hide();
 
$("#eq15text").hide();
 
$("#eq16text").hide();
 
$("#eq17text").hide();
 
$("#eq18text").hide();
 
$("#eq19text").hide();
 
$("#eq20text").hide();
 
$("#eq21text").hide();
 
$("#eq22text").hide();
 
$("#eq23text").hide();
 
$("#eq24text").hide();
 
$("#eq25text").hide();
 
$("#eq26text").hide();
 
$("#eq27text").hide();
 
$("#eq28text").hide();
 
});
 
 
$( "#eq1" ).mouseleave(function() {
 
$(".sidebarright").hide();
 
});
 
 
$( "#eq2" ).mouseenter(function() {
 
$(".sidebarright").show();
 
$("#eq1text").hide();
 
$("#eq2text").show();
 
$("#eq3text").hide();
 
$("#eq4text").hide();
 
$("#eq5text").hide();
 
$("#eq6text").hide();
 
$("#eq7text").hide();
 
$("#eq8text").hide();
 
$("#eq9text").hide();
 
$("#eq10text").hide();
 
$("#eq11text").hide();
 
$("#eq12text").hide();
 
$("#eq13text").hide();
 
$("#eq14text").hide();
 
$("#eq15text").hide();
 
$("#eq16text").hide();
 
$("#eq17text").hide();
 
$("#eq18text").hide();
 
$("#eq19text").hide();
 
$("#eq20text").hide();
 
$("#eq21text").hide();
 
$("#eq22text").hide();
 
$("#eq23text").hide();
 
$("#eq24text").hide();
 
$("#eq25text").hide();
 
$("#eq26text").hide();
 
$("#eq27text").hide();
 
$("#eq28text").hide();
 
});
 
 
$( "#eq2" ).mouseleave(function() {
 
$(".sidebarright").hide();
 
});
 
 
$( "#eq3" ).mouseenter(function() {
 
$(".sidebarright").show();
 
$("#eq1text").hide();
 
$("#eq2text").hide();
 
$("#eq3text").show();
 
$("#eq4text").hide();
 
$("#eq5text").hide();
 
$("#eq6text").hide();
 
$("#eq7text").hide();
 
$("#eq8text").hide();
 
$("#eq9text").hide();
 
$("#eq10text").hide();
 
$("#eq11text").hide();
 
$("#eq12text").hide();
 
$("#eq13text").hide();
 
$("#eq14text").hide();
 
$("#eq15text").hide();
 
$("#eq16text").hide();
 
$("#eq17text").hide();
 
$("#eq18text").hide();
 
$("#eq19text").hide();
 
$("#eq20text").hide();
 
$("#eq21text").hide();
 
$("#eq22text").hide();
 
$("#eq23text").hide();
 
$("#eq24text").hide();
 
$("#eq25text").hide();
 
$("#eq26text").hide();
 
$("#eq27text").hide();
 
$("#eq28text").hide();
 
});
 
 
$( "#eq3" ).mouseleave(function() {
 
$(".sidebarright").hide();
 
});
 
 
$( "#eq4" ).mouseenter(function() {
 
$(".sidebarright").show();
 
$("#eq1text").hide();
 
$("#eq2text").hide();
 
$("#eq3text").hide();
 
$("#eq4text").show();
 
$("#eq5text").hide();
 
$("#eq6text").hide();
 
$("#eq7text").hide();
 
$("#eq8text").hide();
 
$("#eq9text").hide();
 
$("#eq10text").hide();
 
$("#eq11text").hide();
 
$("#eq12text").hide();
 
$("#eq13text").hide();
 
$("#eq14text").hide();
 
$("#eq15text").hide();
 
$("#eq16text").hide();
 
$("#eq17text").hide();
 
$("#eq18text").hide();
 
$("#eq19text").hide();
 
$("#eq20text").hide();
 
$("#eq21text").hide();
 
$("#eq22text").hide();
 
$("#eq23text").hide();
 
$("#eq24text").hide();
 
$("#eq25text").hide();
 
$("#eq26text").hide();
 
$("#eq27text").hide();
 
$("#eq28text").hide();
 
});
 
 
$( "#eq4" ).mouseleave(function() {
 
$(".sidebarright").hide();
 
});
 
 
$( "#eq5" ).mouseenter(function() {
 
$(".sidebarright").show();
 
$("#eq1text").hide();
 
$("#eq2text").hide();
 
$("#eq3text").hide();
 
$("#eq4text").hide();
 
$("#eq5text").show();
 
$("#eq6text").hide();
 
$("#eq7text").hide();
 
$("#eq8text").hide();
 
$("#eq9text").hide();
 
$("#eq10text").hide();
 
$("#eq11text").hide();
 
$("#eq12text").hide();
 
$("#eq13text").hide();
 
$("#eq14text").hide();
 
$("#eq15text").hide();
 
$("#eq16text").hide();
 
$("#eq17text").hide();
 
$("#eq18text").hide();
 
$("#eq19text").hide();
 
$("#eq20text").hide();
 
$("#eq21text").hide();
 
$("#eq22text").hide();
 
$("#eq23text").hide();
 
$("#eq24text").hide();
 
$("#eq25text").hide();
 
$("#eq26text").hide();
 
$("#eq27text").hide();
 
$("#eq28text").hide();
 
});
 
 
$( "#eq5" ).mouseleave(function() {
 
$(".sidebarright").hide();
 
});
 
 
$( "#eq6" ).mouseenter(function() {
 
$(".sidebarright").show();
 
$("#eq1text").hide();
 
$("#eq2text").hide();
 
$("#eq3text").hide();
 
$("#eq4text").hide();
 
$("#eq5text").hide();
 
$("#eq6text").show();
 
$("#eq7text").hide();
 
$("#eq8text").hide();
 
$("#eq9text").hide();
 
$("#eq10text").hide();
 
$("#eq11text").hide();
 
$("#eq12text").hide();
 
$("#eq13text").hide();
 
$("#eq14text").hide();
 
$("#eq15text").hide();
 
$("#eq16text").hide();
 
$("#eq17text").hide();
 
$("#eq18text").hide();
 
$("#eq19text").hide();
 
$("#eq20text").hide();
 
$("#eq21text").hide();
 
$("#eq22text").hide();
 
$("#eq23text").hide();
 
$("#eq24text").hide();
 
$("#eq25text").hide();
 
$("#eq26text").hide();
 
$("#eq27text").hide();
 
$("#eq28text").hide();
 
});
 
 
$( "#eq6" ).mouseleave(function() {
 
$(".sidebarright").hide();
 
});
 
 
$( "#eq7" ).mouseenter(function() {
 
$(".sidebarright").show();
 
$("#eq1text").hide();
 
$("#eq2text").hide();
 
$("#eq3text").hide();
 
$("#eq4text").hide();
 
$("#eq5text").hide();
 
$("#eq6text").hide();
 
$("#eq7text").show();
 
$("#eq8text").hide();
 
$("#eq9text").hide();
 
$("#eq10text").hide();
 
$("#eq11text").hide();
 
$("#eq12text").hide();
 
$("#eq13text").hide();
 
$("#eq14text").hide();
 
$("#eq15text").hide();
 
$("#eq16text").hide();
 
$("#eq17text").hide();
 
$("#eq18text").hide();
 
$("#eq19text").hide();
 
$("#eq20text").hide();
 
$("#eq21text").hide();
 
$("#eq22text").hide();
 
$("#eq23text").hide();
 
$("#eq24text").hide();
 
$("#eq25text").hide();
 
$("#eq26text").hide();
 
$("#eq27text").hide();
 
$("#eq28text").hide();
 
});
 
 
$( "#eq7" ).mouseleave(function() {
 
$(".sidebarright").hide();
 
});
 
 
$( "#eq8" ).mouseenter(function() {
 
$(".sidebarright").show();
 
$("#eq1text").hide();
 
$("#eq2text").hide();
 
$("#eq3text").hide();
 
$("#eq4text").hide();
 
$("#eq5text").hide();
 
$("#eq6text").hide();
 
$("#eq7text").hide();
 
$("#eq8text").show();
 
$("#eq9text").hide();
 
$("#eq10text").hide();
 
$("#eq11text").hide();
 
$("#eq12text").hide();
 
$("#eq13text").hide();
 
$("#eq14text").hide();
 
$("#eq15text").hide();
 
$("#eq16text").hide();
 
$("#eq17text").hide();
 
$("#eq18text").hide();
 
$("#eq19text").hide();
 
$("#eq20text").hide();
 
$("#eq21text").hide();
 
$("#eq22text").hide();
 
$("#eq23text").hide();
 
$("#eq24text").hide();
 
$("#eq25text").hide();
 
$("#eq26text").hide();
 
$("#eq27text").hide();
 
$("#eq28text").hide();
 
});
 
 
$( "#eq8" ).mouseleave(function() {
 
$(".sidebarright").hide();
 
});
 
 
$( "#eq9" ).mouseenter(function() {
 
$(".sidebarright").show();
 
$("#eq1text").hide();
 
$("#eq2text").hide();
 
$("#eq3text").hide();
 
$("#eq4text").hide();
 
$("#eq5text").hide();
 
$("#eq6text").hide();
 
$("#eq7text").hide();
 
$("#eq8text").hide();
 
$("#eq9text").show();
 
$("#eq10text").hide();
 
$("#eq11text").hide();
 
$("#eq12text").hide();
 
$("#eq13text").hide();
 
$("#eq14text").hide();
 
$("#eq15text").hide();
 
$("#eq16text").hide();
 
$("#eq17text").hide();
 
$("#eq18text").hide();
 
$("#eq19text").hide();
 
$("#eq20text").hide();
 
$("#eq21text").hide();
 
$("#eq22text").hide();
 
$("#eq23text").hide();
 
$("#eq24text").hide();
 
$("#eq25text").hide();
 
$("#eq26text").hide();
 
$("#eq27text").hide();
 
$("#eq28text").hide();
 
$(".sidebarleft").show();
 
$("#eq9textleft").show();
 
$("#eq16textleft").hide();
 
$("#eq21textleft").hide();
 
$("#eq22textleft").hide();
 
$("#eq23textleft").hide();
 
});
 
 
$( "#eq9" ).mouseleave(function() {
 
$(".sidebarright").hide();
 
$(".sidebarleft").hide();
 
});
 
 
$( "#eq10" ).mouseenter(function() {
 
$(".sidebarright").show();
 
$("#eq1text").hide();
 
$("#eq2text").hide();
 
$("#eq3text").hide();
 
$("#eq4text").hide();
 
$("#eq5text").hide();
 
$("#eq6text").hide();
 
$("#eq7text").hide();
 
$("#eq8text").hide();
 
$("#eq9text").hide();
 
$("#eq10text").show();
 
$("#eq11text").hide();
 
$("#eq12text").hide();
 
$("#eq13text").hide();
 
$("#eq14text").hide();
 
$("#eq15text").hide();
 
$("#eq16text").hide();
 
$("#eq17text").hide();
 
$("#eq18text").hide();
 
$("#eq19text").hide();
 
$("#eq20text").hide();
 
$("#eq21text").hide();
 
$("#eq22text").hide();
 
$("#eq23text").hide();
 
$("#eq24text").hide();
 
$("#eq25text").hide();
 
$("#eq26text").hide();
 
$("#eq27text").hide();
 
$("#eq28text").hide();
 
});
 
 
$( "#eq10" ).mouseleave(function() {
 
$(".sidebarright").hide();
 
});
 
 
$( "#eq11" ).mouseenter(function() {
 
$(".sidebarright").show();
 
$("#eq1text").hide();
 
$("#eq2text").hide();
 
$("#eq3text").hide();
 
$("#eq4text").hide();
 
$("#eq5text").hide();
 
$("#eq6text").hide();
 
$("#eq7text").hide();
 
$("#eq8text").hide();
 
$("#eq9text").hide();
 
$("#eq10text").hide();
 
$("#eq11text").show();
 
$("#eq12text").hide();
 
$("#eq13text").hide();
 
$("#eq14text").hide();
 
$("#eq15text").hide();
 
$("#eq16text").hide();
 
$("#eq17text").hide();
 
$("#eq18text").hide();
 
$("#eq19text").hide();
 
$("#eq20text").hide();
 
$("#eq21text").hide();
 
$("#eq22text").hide();
 
$("#eq23text").hide();
 
$("#eq24text").hide();
 
$("#eq25text").hide();
 
$("#eq26text").hide();
 
$("#eq27text").hide();
 
$("#eq28text").hide();
 
});
 
 
$( "#eq11" ).mouseleave(function() {
 
$(".sidebarright").hide();
 
});
 
 
$( "#eq12" ).mouseenter(function() {
 
$(".sidebarright").show();
 
$("#eq1text").hide();
 
$("#eq2text").hide();
 
$("#eq3text").hide();
 
$("#eq4text").hide();
 
$("#eq5text").hide();
 
$("#eq6text").hide();
 
$("#eq7text").hide();
 
$("#eq8text").hide();
 
$("#eq9text").hide();
 
$("#eq10text").hide();
 
$("#eq11text").hide();
 
$("#eq12text").show();
 
$("#eq13text").hide();
 
$("#eq14text").hide();
 
$("#eq15text").hide();
 
$("#eq16text").hide();
 
$("#eq17text").hide();
 
$("#eq18text").hide();
 
$("#eq19text").hide();
 
$("#eq20text").hide();
 
$("#eq21text").hide();
 
$("#eq22text").hide();
 
$("#eq23text").hide();
 
$("#eq24text").hide();
 
$("#eq25text").hide();
 
$("#eq26text").hide();
 
$("#eq27text").hide();
 
$("#eq28text").hide();
 
});
 
 
$( "#eq12" ).mouseleave(function() {
 
$(".sidebarright").hide();
 
});
 
 
$( "#eq13" ).mouseenter(function() {
 
$(".sidebarright").show();
 
$("#eq1text").hide();
 
$("#eq2text").hide();
 
$("#eq3text").hide();
 
$("#eq4text").hide();
 
$("#eq5text").hide();
 
$("#eq6text").hide();
 
$("#eq7text").hide();
 
$("#eq8text").hide();
 
$("#eq9text").hide();
 
$("#eq10text").hide();
 
$("#eq11text").hide();
 
$("#eq12text").hide();
 
$("#eq13text").show();
 
$("#eq14text").hide();
 
$("#eq15text").hide();
 
$("#eq16text").hide();
 
$("#eq17text").hide();
 
$("#eq18text").hide();
 
$("#eq19text").hide();
 
$("#eq20text").hide();
 
$("#eq21text").hide();
 
$("#eq22text").hide();
 
$("#eq23text").hide();
 
$("#eq24text").hide();
 
$("#eq25text").hide();
 
$("#eq26text").hide();
 
$("#eq27text").hide();
 
$("#eq28text").hide();
 
});
 
 
$( "#eq13" ).mouseleave(function() {
 
$(".sidebarright").hide();
 
});
 
 
$( "#eq14" ).mouseenter(function() {
 
$(".sidebarright").show();
 
$("#eq1text").hide();
 
$("#eq2text").hide();
 
$("#eq3text").hide();
 
$("#eq4text").hide();
 
$("#eq5text").hide();
 
$("#eq6text").hide();
 
$("#eq7text").hide();
 
$("#eq8text").hide();
 
$("#eq9text").hide();
 
$("#eq10text").hide();
 
$("#eq11text").hide();
 
$("#eq12text").hide();
 
$("#eq13text").hide();
 
$("#eq14text").show();
 
$("#eq15text").hide();
 
$("#eq16text").hide();
 
$("#eq17text").hide();
 
$("#eq18text").hide();
 
$("#eq19text").hide();
 
$("#eq20text").hide();
 
$("#eq21text").hide();
 
$("#eq22text").hide();
 
$("#eq23text").hide();
 
$("#eq24text").hide();
 
$("#eq25text").hide();
 
$("#eq26text").hide();
 
$("#eq27text").hide();
 
$("#eq28text").hide();
 
});
 
 
$( "#eq14" ).mouseleave(function() {
 
$(".sidebarright").hide();
 
});
 
 
$( "#eq15" ).mouseenter(function() {
 
$(".sidebarright").show();
 
$("#eq1text").hide();
 
$("#eq2text").hide();
 
$("#eq3text").hide();
 
$("#eq4text").hide();
 
$("#eq5text").hide();
 
$("#eq6text").hide();
 
$("#eq7text").hide();
 
$("#eq8text").hide();
 
$("#eq9text").hide();
 
$("#eq10text").hide();
 
$("#eq11text").hide();
 
$("#eq12text").hide();
 
$("#eq13text").hide();
 
$("#eq14text").hide();
 
$("#eq15text").show();
 
$("#eq16text").hide();
 
$("#eq17text").hide();
 
$("#eq18text").hide();
 
$("#eq19text").hide();
 
$("#eq20text").hide();
 
$("#eq21text").hide();
 
$("#eq22text").hide();
 
$("#eq23text").hide();
 
$("#eq24text").hide();
 
$("#eq25text").hide();
 
$("#eq26text").hide();
 
$("#eq27text").hide();
 
$("#eq28text").hide();
 
});
 
 
$( "#eq15" ).mouseleave(function() {
 
$(".sidebarright").hide();
 
});
 
 
$( "#eq16" ).mouseenter(function() {
 
$(".sidebarright").show();
 
$("#eq1text").hide();
 
$("#eq2text").hide();
 
$("#eq3text").hide();
 
$("#eq4text").hide();
 
$("#eq5text").hide();
 
$("#eq6text").hide();
 
$("#eq7text").hide();
 
$("#eq8text").hide();
 
$("#eq9text").hide();
 
$("#eq10text").hide();
 
$("#eq11text").hide();
 
$("#eq12text").hide();
 
$("#eq13text").hide();
 
$("#eq14text").hide();
 
$("#eq15text").hide();
 
$("#eq16text").show();
 
$("#eq17text").hide();
 
$("#eq18text").hide();
 
$("#eq19text").hide();
 
$("#eq20text").hide();
 
$("#eq21text").hide();
 
$("#eq22text").hide();
 
$("#eq23text").hide();
 
$("#eq24text").hide();
 
$("#eq25text").hide();
 
$("#eq26text").hide();
 
$("#eq27text").hide();
 
$("#eq28text").hide();
 
$(".sidebarleft").show();
 
$("#eq9textleft").hide();
 
$("#eq16textleft").show();
 
$("#eq21textleft").hide();
 
$("#eq22textleft").hide();
 
$("#eq23textleft").hide();
 
});
 
 
$( "#eq16" ).mouseleave(function() {
 
$(".sidebarright").hide();
 
$(".sidebarleft").hide();
 
});
 
 
$( "#eq17" ).mouseenter(function() {
 
$(".sidebarright").show();
 
$("#eq1text").hide();
 
$("#eq2text").hide();
 
$("#eq3text").hide();
 
$("#eq4text").hide();
 
$("#eq5text").hide();
 
$("#eq6text").hide();
 
$("#eq7text").hide();
 
$("#eq8text").hide();
 
$("#eq9text").hide();
 
$("#eq10text").hide();
 
$("#eq11text").hide();
 
$("#eq12text").hide();
 
$("#eq13text").hide();
 
$("#eq14text").hide();
 
$("#eq15text").hide();
 
$("#eq16text").hide();
 
$("#eq17text").show();
 
$("#eq18text").hide();
 
$("#eq19text").hide();
 
$("#eq20text").hide();
 
$("#eq21text").hide();
 
$("#eq22text").hide();
 
$("#eq23text").hide();
 
$("#eq24text").hide();
 
$("#eq25text").hide();
 
$("#eq26text").hide();
 
$("#eq27text").hide();
 
$("#eq28text").hide();
 
});
 
 
$( "#eq17" ).mouseleave(function() {
 
$(".sidebarright").hide();
 
});
 
 
$( "#eq18" ).mouseenter(function() {
 
$(".sidebarright").show();
 
$("#eq1text").hide();
 
$("#eq2text").hide();
 
$("#eq3text").hide();
 
$("#eq4text").hide();
 
$("#eq5text").hide();
 
$("#eq6text").hide();
 
$("#eq7text").hide();
 
$("#eq8text").hide();
 
$("#eq9text").hide();
 
$("#eq10text").hide();
 
$("#eq11text").hide();
 
$("#eq12text").hide();
 
$("#eq13text").hide();
 
$("#eq14text").hide();
 
$("#eq15text").hide();
 
$("#eq16text").hide();
 
$("#eq17text").hide();
 
$("#eq18text").show();
 
$("#eq19text").hide();
 
$("#eq20text").hide();
 
$("#eq21text").hide();
 
$("#eq22text").hide();
 
$("#eq23text").hide();
 
$("#eq24text").hide();
 
$("#eq25text").hide();
 
$("#eq26text").hide();
 
$("#eq27text").hide();
 
$("#eq28text").hide();
 
});
 
 
$( "#eq18" ).mouseleave(function() {
 
$(".sidebarright").hide();
 
});
 
 
$( "#eq19" ).mouseenter(function() {
 
$(".sidebarright").show();
 
$("#eq1text").hide();
 
$("#eq2text").hide();
 
$("#eq3text").hide();
 
$("#eq4text").hide();
 
$("#eq5text").hide();
 
$("#eq6text").hide();
 
$("#eq7text").hide();
 
$("#eq8text").hide();
 
$("#eq9text").hide();
 
$("#eq10text").hide();
 
$("#eq11text").hide();
 
$("#eq12text").hide();
 
$("#eq13text").hide();
 
$("#eq14text").hide();
 
$("#eq15text").hide();
 
$("#eq16text").hide();
 
$("#eq17text").hide();
 
$("#eq18text").hide();
 
$("#eq19text").show();
 
$("#eq20text").hide();
 
$("#eq21text").hide();
 
$("#eq22text").hide();
 
$("#eq23text").hide();
 
$("#eq24text").hide();
 
$("#eq25text").hide();
 
$("#eq26text").hide();
 
$("#eq27text").hide();
 
$("#eq28text").hide();
 
});
 
 
$( "#eq19" ).mouseleave(function() {
 
$(".sidebarright").hide();
 
});
 
 
$( "#eq20" ).mouseenter(function() {
 
$(".sidebarright").show();
 
$("#eq1text").hide();
 
$("#eq2text").hide();
 
$("#eq3text").hide();
 
$("#eq4text").hide();
 
$("#eq5text").hide();
 
$("#eq6text").hide();
 
$("#eq7text").hide();
 
$("#eq8text").hide();
 
$("#eq9text").hide();
 
$("#eq10text").hide();
 
$("#eq11text").hide();
 
$("#eq12text").hide();
 
$("#eq13text").hide();
 
$("#eq14text").hide();
 
$("#eq15text").hide();
 
$("#eq16text").hide();
 
$("#eq17text").hide();
 
$("#eq18text").hide();
 
$("#eq19text").hide();
 
$("#eq20text").show();
 
$("#eq21text").hide();
 
$("#eq22text").hide();
 
$("#eq23text").hide();
 
$("#eq24text").hide();
 
$("#eq25text").hide();
 
$("#eq26text").hide();
 
$("#eq27text").hide();
 
$("#eq28text").hide();
 
});
 
 
$( "#eq20" ).mouseleave(function() {
 
$(".sidebarright").hide();
 
});
 
 
$( "#eq21" ).mouseenter(function() {
 
$(".sidebarright").show();
 
$("#eq1text").hide();
 
$("#eq2text").hide();
 
$("#eq3text").hide();
 
$("#eq4text").hide();
 
$("#eq5text").hide();
 
$("#eq6text").hide();
 
$("#eq7text").hide();
 
$("#eq8text").hide();
 
$("#eq9text").hide();
 
$("#eq10text").hide();
 
$("#eq11text").hide();
 
$("#eq12text").hide();
 
$("#eq13text").hide();
 
$("#eq14text").hide();
 
$("#eq15text").hide();
 
$("#eq16text").hide();
 
$("#eq17text").hide();
 
$("#eq18text").hide();
 
$("#eq19text").hide();
 
$("#eq20text").hide();
 
$("#eq21text").show();
 
$("#eq22text").hide();
 
$("#eq23text").hide();
 
$("#eq24text").hide();
 
$("#eq25text").hide();
 
$("#eq26text").hide();
 
$("#eq27text").hide();
 
$("#eq28text").hide();
 
$(".sidebarleft").show();
 
$("#eq9textleft").hide();
 
$("#eq16textleft").hide();
 
$("#eq21textleft").show();
 
$("#eq22textleft").hide();
 
$("#eq23textleft").hide();
 
});
 
 
$( "#eq21" ).mouseleave(function() {
 
$(".sidebarright").hide();
 
$(".sidebarleft").hide();
 
});
 
 
$( "#eq22" ).mouseenter(function() {
 
$(".sidebarright").show();
 
$("#eq1text").hide();
 
$("#eq2text").hide();
 
$("#eq3text").hide();
 
$("#eq4text").hide();
 
$("#eq5text").hide();
 
$("#eq6text").hide();
 
$("#eq7text").hide();
 
$("#eq8text").hide();
 
$("#eq9text").hide();
 
$("#eq10text").hide();
 
$("#eq11text").hide();
 
$("#eq12text").hide();
 
$("#eq13text").hide();
 
$("#eq14text").hide();
 
$("#eq15text").hide();
 
$("#eq16text").hide();
 
$("#eq17text").hide();
 
$("#eq18text").hide();
 
$("#eq19text").hide();
 
$("#eq20text").hide();
 
$("#eq21text").hide();
 
$("#eq22text").show();
 
$("#eq23text").hide();
 
$("#eq24text").hide();
 
$("#eq25text").hide();
 
$("#eq26text").hide();
 
$("#eq27text").hide();
 
$("#eq28text").hide();
 
$(".sidebarleft").show();
 
$("#eq9textleft").hide();
 
$("#eq16textleft").hide();
 
$("#eq21textleft").hide();
 
$("#eq22textleft").show();
 
$("#eq23textleft").hide();
 
});
 
 
$( "#eq22" ).mouseleave(function() {
 
$(".sidebarright").hide();
 
$(".sidebarleft").hide();
 
});
 
 
$( "#eq23" ).mouseenter(function() {
 
$(".sidebarright").show();
 
$("#eq1text").hide();
 
$("#eq2text").hide();
 
$("#eq3text").hide();
 
$("#eq4text").hide();
 
$("#eq5text").hide();
 
$("#eq6text").hide();
 
$("#eq7text").hide();
 
$("#eq8text").hide();
 
$("#eq9text").hide();
 
$("#eq10text").hide();
 
$("#eq11text").hide();
 
$("#eq12text").hide();
 
$("#eq13text").hide();
 
$("#eq14text").hide();
 
$("#eq15text").hide();
 
$("#eq16text").hide();
 
$("#eq17text").hide();
 
$("#eq18text").hide();
 
$("#eq19text").hide();
 
$("#eq20text").hide();
 
$("#eq21text").hide();
 
$("#eq22text").hide();
 
$("#eq23text").show();
 
$("#eq24text").hide();
 
$("#eq25text").hi
 

Revision as of 13:45, 19 November 2015

$C$: concentration [nmol/cl]
$\vec{r}$: position vector [µm]
$t$: time [h]
$D$: diffusion constant [µm^2/h]

$C$: concentration [nmol/cl]
$\vec{r}$: position vector [µm]
$t$: time [h]
$\alpha$: production constant [nmol/h]
$\rho_A$: density of type A bacteria [#/cl]

$C$: concentration [nmol/cl]
$\vec{r}$: position vector [µm]
$t$: time [h]
$k$: degradation constant [1/h]

$C$: concentration [nmol/cl]
$\vec{r}$: position vector [µm]
$t$: time [h]
$D$: diffusion constant [µm^2/h]
$\alpha$: production constant [nmol/h]
$\rho_A$: density of type A bacteria [#/cl]
$k$: degradation constant [1/h]

$\vec{r}$: position vector [µm]
$t$: time [h]
$\vec{F}_{applied}$: applied force [mN]
$\gamma$: frictional coefficient [mN*h/µm]

$\vec{r}$: position vector [µm]
$t$: time [h]
$\chi$: chemotactic sensitivity [µm^2*cl/(h*nmol)]
$S$: chemoattractant concentration [nmol/cl]
$\mu$: bacterial diffusion constant [µm^2/h]
$\vec{W}$: Wiener process [h^(1/2)]
$\kappa$: chemotactic sensitivity constant [-]

$n$: bacteria density [#/cl]
$\vec{r}$: position vector [µm]
$t$: time [h]
$\mu$: bacterial diffusion constant [µm^2/h]
$\chi$: chemotactic sensitivity [µm^2*cl/(h*nmol)]
$S$: chemoattractant concentration [nmol/cl]
$\kappa$: chemotactic sensitivity coefficient [-]

$\chi$: chemotactic sensitivity [µm^2*cl/(h*nmol)]
$S$: chemoattractant concentration [nmol/cl]
$\mu$: bacterial diffusion constant [µm^2/h]
$\kappa$: chemotactic sensitivity constant [-]
$S$: chemoattractant concentration [nmol/cl]
$\vec{r}$: position vector [µm]
$t$: time [h]

$E_p$: interaction potential [nJ]
$r_{ij}$: intercellular distance [µm]
$k$: “spring constant” [mN/µm]
$R_{cutoff}$: cutoff intercellular distance [µm]
$R_0$: equilibrium intercellular distance [µm]
$C$: constant to ensure continuity [nJ]

$C_1=-\frac{1}{2}\cdot k_3 \cdot(R_{cutoff}-R_0)^2$
$C_2=-\frac{1}{2}\cdot k_1 \cdot(\frac{R_0}{2}-\frac{k_1+k_2}{k_1}\cdot \frac{R_0}{2})^2 \\ +\frac{1}{2}\cdot k_2 \cdot(\frac{R_0}{2}-R_0)^2 +C_1$

$\vec{F}_i$: cell-cell force acting on cell i [mN]
$E_p$: interaction potential [nJ]
$r_{ij}$: intercellular distance [µm]
$\vec{e}_{ij}$: unit vector from cell i to cell j [-]
$k$: “spring constant” [mN/µm]
$R_{cutoff}$: cutoff intercellular distance [µm]
$R_0$: equilibrium intercellular distance [µm]

$\vec{r}$: position vector [µm]
$t$: time [h]
$\mu$: bacterial diffusion constant [µm^2/h]
$H$: concentration of AHL [nmol/cl]
$\vec{W}$: Wiener process [h^(1/2)]
$\gamma$: frictional coefficient [mN*h/µm]
$E_p$: interaction potential [nJ]
$r_{ij}$: intercellular distance [µm]
$\vec{e}_{ij}$: unit vector from cell i to cell j
$\chi$: chemotactic sensitivity [µm^2*cl/(h*nmol)]
$L$: concentration of leucine [nmol/cl]

$\chi$: chemotactic sensitivity [µm^2*cl/(h*nmol)]
$L$: concentration of leucine [nmol/cl]
$H$: concentration of AHL [nmol/cl]
$\mu$: bacterial diffusion constant [µm^2/h]
$\kappa$: chemotactic sensitivity constant [-]
$\vec{r}$: position vector [µm]
$t$: time [h]

$\mu$: bacterial diffusion constant [µm^2/h]
$H$: concentration of AHL [nmol/cl]
$\vec{r}$: position vector [µm]
$t$: time [h]

$K$: probability density [-]
$x$: stochastic variable [-]

$K_h$: kernel density [#/µm]
$x$: position [µm]
$h$: bandwidth [µm]

$\rho$: agent density [#/cl]
$\beta$: conversion factor [µm/cl]
$N$: number of agents [#]
$K_h$: kernel density [#/µm]
$x$: position [µm]
$x_i$: position of ith agent [µm]

Conversion factor $\beta$ is needed to express the 1-D density [#/µm] in units used in the PDE model [#/cl]. For the 2-D KDE $\beta$ has units [µm^2/cl]. The value of $\beta$ depends on assumptions made about the dimensions of the physical domain which are not modeled. For example, when modeling a 2-D virtual domain this constant will be high if bacteria are densely distributed along the height of the experimental petri dish.

$f$: function, e.g. concentration [nmol/cl]
$\hat{f}$: interpolated function, e.g. [nmol/cl]
$x$: independent variable, e.g. position [µm]
$x_i$: independent variable for ith measurement, e.g. [µm]

$f$: function, e.g. concentration [nmol/cl]
$\hat{f}$: interpolated function, e.g. [nmol/cl]
$x$: independent variable, e.g. position [µm]
$x_i$: independent variable for measurement (i,j), e.g. [µm]
$y$: independent variable, e.g. position [µm]
$y_j$: independent variable for measurement (i,j), e.g. [µm]
$Q_{i,j}$: symbol for grid point corresponding to measurement (i,j)

$f$: function, e.g. concentration [nmol/cl]
$\hat{f}$: interpolated function, e.g. [nmol/cl]
$x$: independent variable, e.g. position [µm]
$x_i$: independent variable for measurement (i,j), e.g. [µm]
$y$: independent variable, e.g. position [µm]
$y_j$: independent variable for measurement (i,j), e.g. [µm]
$Q_{i,j}$: grid point corresponding to measurement (i,j) [-]

$Re$: Reynolds number [-]
$\rho$: mass density [kg/m^3]
$v$: speed [m/s]
$L$: characteristic length [m]
$\eta$: dynamic viscosity [Pa*s]

$F$: force [N]
$m$: mass [kg]
$x$: position [m]
$t$: time [s]
$\eta$: dynamic viscosity [Pa*s]
$a$: radius of particle [m]
$v$: speed [m/s]

Stokes' law states that the force of viscosity for a spherical object moving through a viscous fluid is given by $F=6\pi\eta av$ with $a$ the radius of the sphere. Importantly, this force is proportional to the velocity $v$.

$v$: speed [m/s]
$v_0$: initial speed [m/s]
$\eta$: dynamic viscosity [Pa*s]
$a$: radius of cell [m]
$m$: mass [kg]
$t$: time [s]
$t_0$: start time [s]

To calculate the mass of a bacterium, we assume that the bacterium is spherical with radius $a$ and the density is twice that of water. Therefore, $m=4/3 \cdot \pi \cdot a^3 \cdot (2 \rho)$.

$\mu$: bacterial diffusion constant [µm^2/h]
$\vec{W}$: Wiener process [h^(1/2)]
$\Delta t$: timestep [h]
$\vec{N}$: normal distribution [-]

$\vec{N}$ indicates a vector of which the components are independent random variables that are normally distributed.

$\mu$: bacterial diffusion constant [µm^2/h]
$\Delta t$: timestep [h]
$\vec{N}$: normal distribution [-]

$X_1$: coordinate at next timestep [µm]
$x_1$: coordinate at current timestep [µm]
$\mu$: bacterial diffusion constant [µm^2/h]
$\Delta t$: timestep [h]
$N$: normal distribution [-]

$R$: intercellular distance at next timestep [µm]
$X_1$: coordinate at next timestep [µm]

$X_1$: coordinate at next timestep [µm]
$x_1$: coordinate at current timestep [µm]
$\mu$: bacterial diffusion constant [µm^2/h]
$\Delta t$: timestep [h]
$N$: normal distribution [-]

The Hybrid Model

Introduction

The hybrid model represents an intermediate level of detail in between the colony level model and the internal model. Bacteria are treated as individual agents that behave according to Keller-Segel type stochastic differential equations, while chemical species are modeled using partial differential equations. These different models are implemented and coupled within a single hybrid modeling framework.

Partial Differential Equations

Spatial reaction-diffusion models that rely on Partial Differential Equations (PDEs) are based on the assumption that the collective behavior of individual entities, such as molecules or bacteria, can be abstracted to the behavior of a continuous field that represents the density of those entities. The brownian motion of molecules, for instance, is the result of inherently stochastic processes that take place at the individual molecule level, but is modeled at the density level by Fick’s laws of diffusion. These PDE-based models provide a robust method to predict the evolution of large-scale systems, but are only valid when the spatiotemporal scale is sufficiently large to neglect small-scale stochastic fluctuations and physical granularity. Moreover, such a continuous field approximation can only be made if the behavior of the individual entities is well described.

Agent-Based Models

Agent-based models on the other hand explicitly treat the entities as individual “agents” that behave according to a set of “agent rules”. An agent is an object that acts independently from other agents and is influenced only by its local environment. The goal in agent-based models is to study the emergent systems-level properties of a collection of individual agents that follow relatively simple rules. In smoothed particle hydrodynamics for example, fluids are simulated by calculating the trajectory of each individual fluid particle at every timestep. Fluid properties such as the momentum at a certain point can then be sampled by taking a weighted sum of the momenta of the surrounding fluid particles. A large advantage of agent-based models is that the agent rules are arbitrarily complex and thus they allow us to model systems that do not correspond to any existing or easily derivable PDE model. However, because every agent is stored in memory and needs to be processed individually, simulating agent-based models can be computationally intensive.

Hybrid Modeling Framework

In our system, there are both bacteria and chemical species that spread out and interact on a petri dish to form patterns. On the one hand, the bacteria are rather complex entities that move along chemical gradients and interact with one another. Therefore they are ideally modeled using an agent-based model. On the other hand, the diffusion and dynamics of the chemicals leucine and AHL are easily described by well-established PDEs. To make use of the advantages of each modeling approach, we decided to combine these two different types of modeling in a hybrid modeling framework. In this framework we modeled the bacteria as agents, while the chemical species were modeled using PDEs. There were two challenges to our hybrid approach, namely coupling the models and matching them. Coupling refers to the transfer of information between the models and matching refers to dealing with different spatial and temporal scales to achieve accurate, yet computationally tractable simulations.

In the following paragraphs we first introduce our hybrid model and its coupling. Once the basic framework is established, the agent-based module and PDE module are discussed in more depth and the issue of matching is highlighted. We also expand on important aspects of the model and its implementation such as boundary conditions and choice of timesteps. Then the results for the 1-D model and 2-D model simulations are shown and summarized. Finally, the incorporation of the internal model into the hybrid model is discussed and a proof of concept is demonstrated.

Model Description

System

The main protagonists in our pattern-forming system are cell types A and B, AHL and leucine. Cells A produce AHL as well as leucine. They are unaffected by leucine, while cells B are repelled by leucine. AHL modulates the motility of both cell types A and B, but in opposite ways. High concentrations of AHL will render cell type A unable to swim but will activate cell type B’s motility. Conversely, low concentrations of AHL causes swimming of cell type A and incessant tumbling (thus immobility) of cell type B. Lastly, cells A express the adhesin membrane protein, which causes them to stick to each other. Simply said, our system should produce spots of immobile, sticky groups of A type cells, surrounded by rings of B type cells. Any cell that finds itself outside of the region that it should be in, is able to swim to their correct place and becomes immobile there. More details can be found in the research section.

Partial Differential Equations

As discussed in the previous paragraph, our hybrid model incorporates chemical species using PDEs. In our system these are AHL and leucine. The diffusion of AHL and leucine can be described by the heat equation (1).

$$\frac{\partial C(\vec{r},t)}{\partial t}=D \cdot \nabla^2 C(\vec{r},t) \;\;\; \text{(1)}$$

By using (1) we assume that the diffusion speed is isotropic, i.e. the same in all spatial directions. This also explains why it is called the heat equation, since heat diffuses equally fast in all directions. A detailed explanation of the heat equation can be found in box 1. The second factor that needs be taken into account is the production of AHL and leucine by type A bacteria. In principle, AHL and leucine production is dependent on the dynamically-evolving internal states of all cells of type A. However, for our hybrid model we ignored the inner life of all bacteria and instead assumed that AHL and leucine production is directly proportional to the density of A type cells (2).

$$ \frac{\partial C(\vec{r},t)}{\partial t}=\alpha \cdot \rho_A(\vec{r},t) \;\;\; \text{(2)}$$

In the last paragraph we will reconsider this assumption and assign each cell an internal model. Finally, AHL and leucine are organic molecules and therefore they will degrade over time. We assume first-order kinetics meaning that the rate at which AHL and similarly leucine disappear is proportional to their respective concentrations (3a and 3b) assuming neutral pH [6].

$$ \frac{\partial C_{AHL}(\vec{r},t)}{\partial t}=-k_{AHL}\cdot C_{AHL}(\vec{r},t) \;\;\; \text{(3a)} $$ $$ \frac{\partial C_{leucine}(\vec{r},t)}{\partial t}=-k_{leucine}\cdot C_{leucine}(\vec{r},t) \;\;\; \text{(3b)} $$

Putting it all together, we obtain (4), both for AHL and leucine.

$$ \frac{\partial C(\vec{r},t)}{\partial t}=D \cdot \nabla^2 C(\vec{r},t)+\alpha \cdot \rho_A(\vec{r},t)-k\cdot C(\vec{r},t) \;\;\; \text{(4)} $$

Note that these equations have exactly the same form as the equations for AHL and leucine in the colony level model. The crucial difference however lies in the calculation of the density of cells of type A. In contrast to the colony level model, in this model the cell density is not calculated explicitly with a PDE and is therefore not trivially known. Therefore a method to extract a density field from a spatial distribution of agents is necessary. This is addressed in the subparagraph below on coupling.

Box 1: Heat Equation

Heat Equation

The left-hand side of the heat equation (1) represents the rate of accumulation of a chemical, while the right-hand side is proportional to the Laplacian of its concentration field, which is a second-order differential operator. This equation can be easily understood by considering a one-dimensional concentration profile, as shown in Figure 1: if the concentration can be approximated as a convex parabolic function, the second derivative is positive and therefore the rate of accumulation is positive (i.e. more accumulation). If on the other hand the concentration resembles a concave parabolic function, the second derivative is negative and the rate of accumulation as well (i.e. depletion). A special case occurs when the concentration profile takes on a linear form. Everything that moves into the point goes out at the other side and as a result there is no net accumulation over time.

Illustration of the heat equation

Figure 1
Illustration of the heat equation. Click to enlarge.

In the videoplayer below we demonstrate how the heat equation smoothes out an initially heterogeneous concentration profile. When only diffusion is acting on the system, it will always evolve to a uniformly flat concentration profile, regardless of the initial conditions.


Agents

To model bacteria movement on the other hand, we used an agent-based model that explicitly stored individual bacteria as agents. Spatial coordinates are associated with each agent, specifying their location. After solving the equation of motion for all agents based on their environment, these coordinates are updated at every timestep. In principle, Newton’s second law of motion has to be solved for all bacteria. However, since bacteria live in a low Reynolds (high friction) environment, the inertia of the bacteria can be neglected. This is because an applied force will immediately be balanced out by an opposing frictional force, with no noticeable acceleration or deceleration phase taking place. This eliminates the inertial term and simplifies Newton’s second law to (5).

$$ \frac{d^2 \vec{r}(t)}{dt^2}=\sum_{i} \vec{F}_{applied,i}-\gamma \cdot \frac{d \vec{r}(t)}{dt}=0 $$ $$\Rightarrow \frac{d \vec{r}(t)}{dt}=\frac{1}{\gamma} \cdot \sum_{i} \vec{F}_{applied,i} \;\;\; \text{(5)} $$

Basically, the velocity can be calculated as the sum of all applied forces times divided by a frictional coefficient. For more info about the Reynolds number and “life at low Reynolds number”, we refer to box 2. In the following paragraphs we will investigate the different forces acting on the bacteria and ultimately superimpose them to obtain the final equation of motion.

Box 2: Life at Low Reynolds Number

Life at Low Reynolds Number

The Reynolds number in fluid mechanics is a dimensionless number that characterizes different flow regimes. It is most commonly used to determine whether laminar or turbulent flow will take place in a hydrodynamic system (Figure 2).

Flow regimes for different Reynolds numbers

Figure 2
Flow regimes for different Reynolds numbers. Click to enlarge.

In general however, it quantifies the ratio of inertial forces to viscous forces and is defined as (B2.1).

$$ Re=\frac{\text{inertial forces}}{\text{viscous forces}} =\frac{\rho v^2 L^2}{\eta v L} =\frac{\rho v L}{\eta} \;\;\; \text{(B2.1)} $$

For example, if the Reynolds number is high, the inertia of fluids in motion dominate and turbulent flow will occur. When it is low however, the viscous forces dampen the kinetic energy of fluid particles and stabilize the flow profile, ultimately achieving a regular, laminar flow.

The reason we mention the Reynolds number however is not to study the flow of fluids, but to characterize the behavior of swimming objects inside of a stationary fluid. When applied to objects in a fluid, the Reynolds number tells us whether their inertia can be neglected or not. For bacteria, the characteristic length $L$ is on the order of micrometers, which is quite small. Therefore the Reynolds number is always very small and hence viscous forces dominate the motion of bacteria. This is what allows us to eliminate the inertial term in Newton's second law and greatly simplify the equation of motion.

To justify this simplification we will go through a numerical example. Take a bacterium of size $L = 2 a = 1 \cdot \mu m$, swimming through water with a speed of $v=20 \cdot \mu m/s$. Water has a density of $\rho = 1 \; 000 \cdot kg/m^3$ and a viscosity of $\eta = 0.001 \cdot Pa/s$.

Sketch of swimming bacterium

Figure 3
Sketch of swimming bacterium. Click to enlarge.

Putting these values in (B2.1) yields a Reynolds number of $Re = 2 \cdot 10^{-5} << 1$. Clearly we are in an extremely low Reynolds number regime. To show what this really means, suppose the bacterium stops propelling itself. How long will it continue to move relying only on its inertia? Assuming Stokes' law, we obtain (B2.2) as the equation of motion.

$$ F=m\cdot \frac{d^2x}{dt^2} \;\;\; \text{(B2.2a)} $$ $$ -6 \pi \eta a v=m\cdot \frac{dv}{dt} \;\;\; \text{(B2.2b)} $$

Solving this differential equation yields (B2.3).

$$ v=v_0 \cdot \text{exp} \Bigg[ -\frac{6 \pi \eta a}{m} \cdot (t-t_0) \Bigg] \;\;\; \text{(B2.3)} $$

The characteristic time constant is $\tau = 6\pi\eta a/m \approx 0.1 \cdot \mu s$, from which we calculate that the total distance traversed is $\Delta x < 2 \cdot pm$. This coasting distance is 6 orders of magnitude smaller than its size. Moreover, the time spent coasting is extremely short. Thus, once the bacterium stops propelling itself, it is safe to assume that whatever kinetic energy it had is immediately absorbed by hydrodynamic friction, instantly halting the bacterium. Therefore, we can neglect the inertial term in Newton's second law (5).

Stochastic Differential Equation

When we try to calculate the physical “chemotactic force”, propelling bacteria towards chemoattractants or away from chemorepellents, we find that it is not easily measured nor derived. Therefore, as a workaround we base the equation of motion on (6), a stochastic differential equation (SDE) that describes the motion of a single particle in a N-particle system that is governed by a Keller-Segel type PDE in the limit of $N \rightarrow \infty$. This Keller-Segel type PDE (7) describes the evolution of a bacteria density field $n$ moving towards some chemoattractant $S$.

$$ d\vec{r}_i(t)=\chi (S) \cdot \nabla S(\vec{r},t)\cdot dt + \sqrt{2 \cdot \mu}\cdot d\vec{W} \;\;\; \text{(6a)} $$ $$ d\vec{r}_i(t)=\mu \cdot \frac{\kappa}{S(\vec{r},t)} \cdot \nabla S(\vec{r},t)\cdot dt + \sqrt{2 \cdot \mu}\cdot d\vec{W} \;\;\; \text{(6b)} $$

$$ \frac{\partial n(\vec{r},t)}{\partial t}=\mu \cdot \nabla^2 n - \nabla (n \cdot \chi(S) \cdot\nabla S(\vec{r},t)) \;\;\; \text{(7a)} $$ $$ \frac{\partial n(\vec{r},t)}{\partial t}=\mu \cdot \nabla^2 n - \nabla (n \cdot \mu \cdot \frac{\kappa}{S(\vec{r},t)} \cdot\nabla S(\vec{r},t)) \;\;\; \text{(7b)} $$

Put differently, when infinitely many particles obey (6), they will exhibit Keller-Segel type spatial dynamics as in (7). In a sense, we’re using a “reverse-engineered” particle equation that corresponds to the Keller-Segel field equation. A detailed theoretical treatment of (6) is outside the scope of this model description because it contains a stochastic variable. The traditional rules of calculus do not apply anymore for stochastic differential equations and a different mathematical theory called Itô calculus is required. It is sufficient to say that the second term containing $dW$ accounts for Brownian motion in the form of random noise added to the displacement of the agents, causing them to diffuse, and that it is governed by the diffusion coefficient $\mu$.

The first term in (6) on the other hand is easily understood as an advective or drift term (net motion) dependent on S, pushing the agents along a positive gradient (for negative chemotaxis the sign is reversed). The chemotactic drift hence points towards an increasing concentration of the chemoattractant. The advective properties are governed by the chemotactic sensitivity function $\chi (S)$. For our model we used a function of the form (8).

$$ \chi(S)=\mu \cdot \frac{\kappa}{S(\vec{r},t)} \;\;\; \text{(8)} $$

The first important thing to note is that we assume $\chi (S)$ to be proportional to $1/S$. This is because Keller and Segel proved that their corresponding PDE model only yields travelling wave solutions if $\chi (S)$ contains such a singularity at some critical concentration $S_{crit}$, and multiplying by $1/S$ is the simplest way to introduce a singularity at $S_{crit} = 0$. Secondly, the proportionality constant is composed of two factors, namely the bacterial diffusion coefficient $\mu$ and chemotactic sensitivity constant $\kappa$. This is done for two reasons. Firstly, when $\mu$ is lowered, both chemotactic and random motion is reduced, which emulates the state of inactivated motility due to high or low concentrations of AHL. Secondly, defining a separate chemotactic sensitivity constant allows us to examine the effect of the relative strength of chemotaxis versus random motion on pattern formation. Conceptually, this term can be viewed as a sort of chemotactic force, defined with regards to a mobility coefficient $\mu$ instead of a frictional coeffient, which are the inverse of each other.

Cell-Cell Interactions

In addition to chemotaxis and diffusion, cell-cell interactions play an important role in pattern formation and also need to be modeled. Bacteria have finite size and therefore multiple bacteria cannot occupy the same space. Moreover, an important mechanism in our system is the aggregation of cells A due to the sticky adhesin protein membrane. To take these mechanisms into account we modeled two types of cell-cell interactions: the purely repulsive interaction of cell B with another cell B and with cell A, and the attractive-repulsive interaction of cell A with another cell A. The interaction between two cells is usually expressed by a potential energy curve defined over the distance between the centers of mass of the two cells (Figure 2).

Illustration of cell-cell interaction potential

Figure 2
Cell-cell interaction potential. Click to enlarge.

Note that the potential energy remains constant after a certain distance, which means that the cells stop interacting. Also, as two cells move closer together, they hit a wall where the potential energy curve abruptly goes to infinity. The reason for this is that two cells cannot occupy the same space and therefore smaller intercellular distances are not allowed. Implementing this theoretical potential is however not possible because the displacement of bacteria is a stochastic process and the bacteria could randomly jump beyond the potential wall, where the force is ill defined. Therefore, we’ve decided to define a piecewise quadratic potential (9), which results in a piecewise linear force that resembles Hooke’s law, but with three different “spring constants” acting in different intervals of intercellular distances (10), as illustrated in Figure 3. The force is defined with respect to the unit vector pointing towards the other cell, meaning that a positive force corresponds to an attractive force and vice-versa.

Cell-cell interaction potential and force curves

Figure 3
Cell-cell interaction potential and force curves. Click to enlarge.

$$ E_{p,attr-rep}(r_{ij})=\left\{\begin{matrix} 0 & R_{cutoff}\leq r_{ij}\\ \frac{1}{2}\cdot k_3 \cdot(r_{ij}-R_0)^2+C_1 & R_0 \leq r_{ij} < R_{cutoff} \\ \frac{1}{2}\cdot k_2 \cdot(r_{ij}-R_0)^2+C_1 & \frac{R_0}{2} \leq r_{ij} < R_0\\ \frac{1}{2}\cdot k_1 \cdot(r_{ij}-\frac{k_1+k_2}{k_1}\cdot \frac{R_0}{2})^2+C_2 & 0 \leq r_{ij} < \frac{R_0}{2} \end{matrix}\right. \;\;\; \text{(9a)}$$ $$ E_{p,rep}(r_{ij})=\left\{\begin{matrix} 0 & R_0\leq r_{ij}\\ \frac{1}{2}\cdot k_2 \cdot(r_{ij}-R_0)^2 & \frac{R_0}{2} \leq r_{ij} < R_0\\ \frac{1}{2}\cdot k_1 \cdot(r_{ij}-\frac{k_1+k_2}{k_1}\cdot \frac{R_0}{2})^2 & 0 \leq r_{ij} < \frac{R_0}{2} \end{matrix}\right. \;\;\; \text{(9b)}$$

$$ \vec{F}_{i,attr-rep}(r_{ij})= \frac{\partial E_{p,attr-rep}(r_{ij})}{\partial r_{ij}} \cdot \vec{e}_{ij}= \left\{\begin{matrix} \vec{0} & R_{cutoff}\leq r_{ij}\\ k_3 \cdot(r_{ij}-R_0) \cdot \vec{e}_{ij} & R_0 \leq r_{ij} < R_{cutoff} \\ k_2 \cdot(r_{ij}-R_0) \cdot \vec{e}_{ij} &\frac{R_0}{2} \leq r_{ij} < R_0\\ k_1 \cdot(r_{ij}-\frac{k_1+k_2}{k_1}\cdot \frac{R_0}{2}) \cdot \vec{e}_{ij} & 0 \leq r_{ij} < \frac{R_0}{2} \end{matrix}\right. \;\;\; \text{(10a)} $$ $$ \vec{F}_{i,rep}(r_{ij})= \frac{\partial E_{p,rep}(r_{ij})}{\partial r_{ij}} \cdot \vec{e}_{ij}= \left\{\begin{matrix} \vec{0} & R_0\leq r_{ij}\\ k_2 \cdot(r_{ij}-R_0) \cdot \vec{e}_{ij} & \frac{R_0}{2} \leq r_{ij} < R_0\\ k_1 \cdot(r_{ij}-\frac{k_1+k_2}{k_1}\cdot \frac{R_0}{2}) \cdot \vec{e}_{ij} & 0 \leq r_{ij} < \frac{R_0}{2} \end{matrix}\right. \;\;\; \text{(10b)} $$

Between A type cells, there is a region of attraction $(R_0 \leq r_{ij} < R_{cutoff})$, where the force points towards the other cell, hence moving them closer together. In the repulsive domain $(r_{ij} < R_0)$, two regions were defined, emulating lower repulsive forces $(\frac{R_0}{2} \leq r_{ij} < R_0)$ and higher repulsive forces due to a higher spring constant when the cells are even closer $(r_{ij} < \frac{R_0}{2})$. For the purely repulsive interaction scheme there is no attraction and therefore the spring constant for $R_0 \leq r_{ij}$ is zero. More details about the implementation of the cell-cell interaction scheme, more specifically regarding the nearest-neighbor search algorithm, can be found in the paragraph on the agent-based module below.

Equation of Motion

Now we are ready to construct the equation of motion for cell type A and B as a superposition of an adapted Keller-Segel SDE (6) and cell-cell interaction forces (10), yielding (11).

$$ d\vec{r}_{A_i}(t)= \sqrt{2 \cdot \mu_A (H)}\cdot d\vec{W} + \frac{1}{\gamma}\cdot\Bigg[ \sum^{A \backslash \{ A_i\}}_j \frac{dE_{p,attr-rep}(r_{ij})}{dr_{ij}} \cdot \vec{e}_{ij} +\sum^{B}_j\frac{dE_{p,rep}(r_{ij})}{dr_{ij}}\cdot \vec{e}_{ij} \Bigg]\cdot dt \;\;\; \text{(11a)} $$ $$ d\vec{r}_{B_i}(t)= \sqrt{2 \cdot \mu_B(H)}\cdot d\vec{W} + \chi(L,H) \cdot \nabla L(\vec{r},t)\cdot dt + \frac{1}{\gamma}\cdot\Bigg[ \sum^{A\cup B\backslash \{ B_i\}}_j \frac{dE_{p,rep}(r_{ij})}{dr_{ij}}\cdot \vec{e}_{ij} \Bigg]\cdot dt \;\;\; \text{(11b)} $$

$$ \chi(L,H)= -\mu_{B}(H) \cdot \frac{\kappa}{L(\vec{r},t)} \;\;\; \text{(12)} $$

$$ \mu_A(H)=\left\{\begin{matrix} \mu_{A,high} & H(\vec{r},t) < H_{A,threshold}\\ \mu_{A,low} & H(\vec{r},t) \geq H_{A,threshold} \end{matrix}\right. \;\;\; \text{(13a)} $$ $$ \mu_B(H)=\left\{\begin{matrix} \mu_{B,high} & H(\vec{r},t) < H_{B,threshold}\\ \mu_{B,low} & H(\vec{r},t) \geq H_{B,threshold} \end{matrix}\right. \;\;\; \text{(13b)} $$

Bacteria of type A are not attracted nor repelled by leucine, so the chemotactic term falls away. All cell-cell forces are summed up to find a net force, taking into account the two different potentials due to the different interaction types. As discussed before, this net force times a constant yields the velocity due to that force, which is then multiplied by $dt$ to obtain the displacement. For B type cells, the chemotactic term models the repulsive chemotaxis away from leucine. The chemotactic sensitivity function (12) has a negative sign signifying that B type cells are repelled by leucine. The cell-cell interaction term in this case is simpler because B type cells only interact repulsively.

Note that the diffusion coefficient of cell types A and B (13) switches based on the local concentration of AHL relative to a threshold AHL value, which simulates the dependency of cellular motility on AHL. For B type cells the cellular motility depends explicitly on AHL due to the synthetic genetic circuit we have built into them. On the other hand, in our model the motility of A type cells should not depend on AHL directly, but since high concentrations of AHL are caused by high densities of aggregated A type cells, the bacteria typically will not be motile in high AHL concentrations because they stick to neighboring cells.

The agent-based module is now fully defined but one crucial issue was skipped: AHL and leucine concentrations are calculated using PDEs and are therefore only known at grid points. Agents on the other hand reside in the space between grid points and require local concentrations as inputs to calculate their next step. This problem is part of the coupling aspect in our hybrid modeling framework and is discussed below.

Coupling

At this point, both the PDE module and agent-based module have been established, but the issue remains that the individual modules take inputs and return outputs which are not directly compatible. In the framework of PDEs, entities are described by density/concentration fields. In the agent-based module however, entities are represented by discrete objects with exact spatial locations. The key challenge posed by the hybrid model is to integrate these different approaches by bidirectionally converting and exchanging information between the modules. Coupling the modules is the essential component that makes the hybrid model hybrid. It allows us to leverage the advantages of both modeling approaches, while circumventing their drawbacks. But although hybrid modeling opens up many new avenues for novel modeling methods, it comes with its own diverse set of issues and peculiarities that need to be addressed before it can be successfully applied.

In the following paragraphs the basic scheme for coupling the PDE and agent-based modules in our model is introduced, after which the theoretical treatment of our hybrid model is complete. The difficulties that arise from linking the modules together are discussed further below in the subsection on matching.

Agent-Based to PDE

As described above, the agents’ effect in the PDE is modeled as a source term that is proportional to the agent density. This approach is essentially the same approach taken in the colony level model for the bacterial production of AHL and leucine. However, in the colony level model the bacteria density is explicitly calculated at the grid points, while the agent-based model essentially considers a set of points scattered in space. A simple first-order approach would be to determine the closest grid point to any agent and simply increment a counter belonging to that grid point. This results in a histogram, which can be used directly to represent the agent density. However, the resulting density is a blocky, nonsmooth function which poorly represents the underlying agent distribution. The effect of a single agent is artificially confined to a single grid point, while in reality an agent’s influence could reach much further than a single grid point. The shape of a histogram is also very dependent on the bin size, which in this case corresponds to the grid spacing so it cannot be independently tuned. To decouple grid spacing and agent density and achieve a smoother density function, we made use of a more sophisticated technique called kernel density estimation (KDE).

KDE is used in statistics to estimate the probability density of a set discrete data derived from a random process. The basic idea consists of defining a kernel function that represents the density of a single data point, then centering kernel functions on every data point and summing them all up to achieve a smooth overall density function, as demonstrated in the Figure 4.

Kernel sum

Figure 4
Kernel density estimation. Click to enlarge.

This kernel function can be anything as long is it continuous, symmetric and integrates to 1, since it represents one data point or one agent in our case. Some of the most common kernel functions include Gaussian kernels, triangular kernels and Epanechnikov kernels (14), which are shown in Figure 5.

Gaussian, triangular and Epanechnikov kernel functions

Figure 5
Gaussian, triangular and Epanechnikov kernel functions. Click to enlarge.

$$ K_{Gaus}(x)=\frac{1}{\sqrt{2 \cdot \pi}} \cdot e^{-\frac{1}{2}x^2} \;\;\; \text{(14a)} $$ $$ K_{tri}(x)=\left\{\begin{matrix} 0 & 1 < |x| \\ 1-|x| & |x| \leq 1 \end{matrix}\right. \;\;\; \text{(14b)} $$ $$ K_{Epa}(x)=\left\{\begin{matrix} 0 & 1 < |x| \\ \frac{3}{4} \cdot (1-x^2) & |x| \leq 1 \end{matrix}\right. \;\;\; \text{(14c)} $$

In practice, scaled versions of standard kernel functions are used, which are of the form (15).

$$ K_h(x)=\frac{1}{h} \cdot K \bigg( \frac{x}{h} \bigg) \;\;\; \text{(15)} $$

Importantly, the scaled functions inherit the kernel function properties, but are either broader or narrower. The degree to which the shape of a kernel function is stretched or squeezed depends on the scaling factor h, which is why it is called the bandwidth, see Figure 6.

Epanechnikov kernel using different bandwidths

Figure 6
Epanechnikov kernels using different bandwidths. Click to enlarge.

This parameter gives us the freedom to define how far the influence of an agent reaches and how smooth the resulting density function looks like. An example is given in Figure 7, where 100 agents were sampled from a normal distribution. The red dotted line indicates the true underlying distribution, while the blue line is the kernel density estimation calculated using different bandwidths. The left graph shows a undersmoothed (small bandwidth) KDE, which oscillates wildly. The right graph on the other hand is oversmoothed (large bandwidth), leading to a misleadingly wide KDE.

Kernel density estimation using different bandwidths

Figure 7
Kernel density estimation using different bandwidths. Click to enlarge.

In summary, using kernel density estimation we can express the agent density in the form of (16), providing the appropriate input for the PDE model.

$$ \rho(x)=\beta \cdot \sum^N_{i} K_h(x-x_i) \;\;\; \text{(16)} $$

The extension to higher dimensions is called multivariate kernel density estimation and is rather non-trivial since the bandwidth is then defined as a matrix instead of a scalar, which allows for many more variations on the same basis kernel function. A detailed analysis of multivariate kernel density estimation will not be given here. It suffices to say that for our 2-D hybrid model we used the simplest possible bandwidth matrix, which is the identity matrix times a constant.

PDE to Agent-Based

The final component of our hybrid model is the mapping of the PDE model to the agent-based model. The latter model works with discrete objects that have continuous coordinates, meaning that they can be located at any point of the domain. As we have seen in (11), the agents need the local concentration of AHL and leucine, as well as the gradient of leucine in order to update their positions. In the PDE model however, the domain is discretized into a grid and concentrations are only defined at grid points. Therefore, in order to transfer information from the PDE model to the agent-based model we need to translate these grid values into values for any given position within the domain. We achieved this for the 1-D hybrid model by taking the two nearest grid points to any agent and employing linear interpolation to derive an approximate local field value. Similarly, for the 2-D hybrid model we took the four nearest grid points and employed bilinear interpolation, as explained in box 3. With the end result (B3.3), the concentration or gradient of AHL or leucine can be derived at any point of the domain. Therefore, the agent-based module has access to the information that it needs from the PDE module.

Box 3: Bilinear Interpolation

Bilinear Interpolation

Interpolation is used to obtain a continuous function when a finite amount of data points are known. Piecewise linear interpolation is the most simple scheme to interpolate 1-dimensional data.

Illustration of linear interpolation

Figure 8
Illustration of linear interpolation. Click to enlarge.

The technique involves defining linear functions connecting consecutive data points to fill in the gaps, using (B3.1).

$$ \hat{f}(x)=f(x_i)+\frac{f(x_{i+1})-f(x_i)}{x_{i+1}-x_i} \cdot (x-x_i) \; \text{for} \; x_i < x < x_{i+1} \;\;\; \text{(B3.1)} $$

Bilinear interpolation is the generalization of this interpolation technique to two dimensions, illustrated using Figure 9.

Illustration of bilinear interpolation

Figure 9
Illustration of bilinear interpolation. Click to enlarge.

Let’s assume we want to interpolate a two-dimensional function within a rectangular domain where the function values are only known on the four adjacent grid points $Q$. First we linearly interpolate in the x-direction on the lines connecting $Q_{i,j}$ and $Q_{i+1,j}$, $Q_{i,j+1}$ and $Q_{i+1,j+1}$, which yields (B3.2).

$$ \hat{f}(x,y_j)=\frac{x_{i+1}-x}{x_{i+1}-x_i} \cdot f(Q_{i,j}) + \frac{x-x_i}{x_{i+1}-x_i} \cdot f(Q_{i+1,j}) \;\;\; \text{(B3.2a)} $$ $$ \hat{f}(x,y_{j+1})=\frac{x_{i+1}-x}{x_{i+1}-x_i} \cdot f(Q_{i,j+1}) + \frac{x-x_i}{x_{i+1}-x_i} \cdot f(Q_{i+1,j+1}) \;\;\; \text{(B3.2b)} $$

We then obtain interpolated values for all positions on the bottom and top edge. The key idea then is to take these values and interpolate in the y-direction for any given x-position, yielding (B3.3).

$$ \hat{f}(x,y)=\frac{y_{j+1}-y}{y_{j+1}-y_j} \cdot \hat{f}(x,y_j) + \frac{y-y_j}{y_{j+1}-y_j} \cdot \hat{f}(x,y_{j+1}) \\ =\frac{y_{j+1}-y}{y_{j+1}-y_j} \cdot \Bigg[ \frac{x_{i+1}-x}{x_{i+1}-x_i} \cdot f(Q_{i,j}) + \frac{x-x_i}{x_{i+1}-x_i} \cdot f(Q_{i+1,j}) \Bigg] \\ +\frac{y-y_j}{y_{j+1}-y_j} \cdot \Bigg[ \frac{x_{i+1}-x}{x_{i+1}-x_i} \cdot f(Q_{i,j+1}) + \frac{x-x_i}{x_{i+1}-x_i} \cdot f(Q_{i+1,j+1}) \Bigg] \;\;\; \text{(B3.3)} $$

This then gives interpolated values for any point within the domain. Note that the resulting function is not linear, but quadratic since it contains the term $x \cdot y$.

Synthesis

The diffusion, production and degradation of AHL and leucine are described by the PDEs (4). At the same time, the random movement, chemotaxis and intercellular interactions of cell types A and B are captured by the stochastic equation of motion (11). Translating the spatial distribution of cells type A into a density field to feed into the PDE module can be accomplished by using kernel density estimation (16). Finally, the agents can request concentrations and gradients at arbitrary positions from the PDE module using bilinear interpolation (B3.3). Taken together, these equations describe the individual modules, as well as the two-way information transfer in between them, and thus they fully define our hybrid model. A graphical summary is shown in Figure 10.

Graphical summary of hybrid model

Figure 10
Graphical summary of hybrid model. Click to enlarge.

Implementation

Agent-Based Module

Introductory text about the Agent-Based Module

Nearest-Neighbor Algorithm

The cell-cell interactions have already been fully described in the paragraphs above. However, solving the equation of motion of an agent in its current form requires the computation of the force due to every other agent. If we take N to be the number of agents, that means $N*(N-1)$ amount of force computations are needed in total and therefore the computation time grows as $O(N^2)$. To put this in perspective, if we simulate a thousand agents, the amount of interactions adds up to around one million. This puts a heavy computational load on the model which can be reduced greatly by using a smarter algorithm. The crucial observation here is that the force goes to zero when the distance between two cells is larger than $2\cdot r_{cutoff}$. An agent therefore only interacts with its nearest neighbors, meaning that the naive implementation wastes a lot of time computing interactions which have no effect anyways. Moreover, as the simulation progresses, an agent’s neighbors do not vary greatly from one timestep to another. A more efficient approach then consists of periodically searching for the nearest neighbors within a fixed radius, storing them in a list for every agent and updating their coordinates for all intermediate timesteps. Since the agents are programmed to repel each other when they approach each other too closely, they will evolve to a rather uniform distribution. The most appropriate fixed-radius nearest neighbor search algorithm in that case is the cell technique. In this algorithm, the agents are structured by placing a rectangular grid of cells over the domain and assigning every agent to a cell (fig9). Determining which cell an agent belongs to is easily done by rounding down the x and y-coordinates to the nearest multiple of the grid spacing. If the grid spacing is chosen so that interactions do not reach further than adjacent cells, every agent only needs to look for neighbors in 9 cells (its own cells plus 8 adjacent cells) instead of the entire domain.

Choice of Timestep

In contrast to PDE models, there is no danger of instability due to a large timestep in agent-based models. However, since the time step determines how far a cell can “jump” away from its previous position, we need to ensure that it is not so large as to negate the effect of intercellular interactions. To understand this, we first discretize the Brownian motion term in the agent’s equation of motion (17).

$$ \sqrt{2 \cdot \mu}\cdot d\vec{W} \rightarrow \sqrt{2 \cdot \mu \cdot \Delta t}\cdot \vec{N}(0,1) \;\;\; \text{(17)} $$

The details of this discretization scheme goes beyond the scope of this text and will not be given here. It can be proven that this term follows a distribution as (18).

$$ \sqrt{2 \cdot \mu \cdot \Delta t}\cdot \vec{N}(0,1) \sim \vec{N}(0,2 \cdot \mu \cdot \Delta t) \;\;\; \text{(18)} $$

Now assume that there are two A type cells who are exactly $2 \cdot r_0$ apart. If the timestep is too large, they might “jump” out of each other’s influence at the next iteration. This situation is highly unrealistic because the cells actually should feel each other’s pull while they randomly drift away from each other. If the timestep is too large, there’s a high probability that the cells detach without experiencing the potential at all and therefore the probability of detachment is artificially inflated. Of course, we can never completely eliminate the possibility of a sudden detachment because the random is sampled by a normal distribution, but we can define a desired probability of sudden detachment $P^-$. Then, for a given $r_0$, $r_{cutoff}$ and $\mu_A$, we can determine how small the timestep has to be for the sudden detachment to happen only with probability $P^-$. First, we assume two A type cells with coordinates $(x_1,y_1)$ and $(x_2,y_2)$ that are at a distance $2 \cdot r_0 = R_0$ from each other. We then define the coordinates at the next timestep as stochastic variables $X_1$, $Y_1$, $X_2$ and $Y_2$, which are distributed as (19).

$$ X_1 \sim x_1 + \sqrt{2 \cdot \mu \cdot \Delta t}\cdot {N}(0,1) \sim {N}(x_1,2 \cdot \mu \cdot \Delta t) \;\;\; \text{(19a)} $$ $$ Y_1 \sim y_1 + \sqrt{2 \cdot \mu \cdot \Delta t}\cdot {N}(0,1) \sim {N}(y_1,2 \cdot \mu \cdot \Delta t) \;\;\; \text{(19b)} $$ $$ X_2 \sim x_2 + \sqrt{2 \cdot \mu \cdot \Delta t}\cdot {N}(0,1) \sim {N}(x_2,2 \cdot \mu \cdot \Delta t) \;\;\; \text{(19c)} $$ $$ Y_2 \sim y_2 + \sqrt{2 \cdot \mu \cdot \Delta t}\cdot {N}(0,1) \sim {N}(y_2,2 \cdot \mu \cdot \Delta t) \;\;\; \text{(19d)} $$

We are not interested in the future coordinates as such, but rather in the future intercellular distance, which is composed of these future coordinates and therefore also is a stochastic variable. For ease of calculation, we work with the square of this variable, which we name as R² and define as (20).

$$ R^2=(X_1-X_2)^2 + (Y_1-Y_2)^2 \;\;\; \text{(20)} $$

It is easy to see that $X_1-X_2$ and $Y_1-Y_2$ are both noncentral (mean is not equal to zero) normally distributed variables as well. Therefore $R^2$ obeys a noncentral chi-squared distribution with 2 degrees of freedom. Before proceeding further, we have to make our $R^2$ conform to a standard noncentral chi-squared distribution. First we normalize $X_1-X_2$ and $Y_1-Y_2$. We observe that $X_1-X_2$ and $Y_1-Y_2$ are distributed as (21).

$$ \begin{split} X_1-X_2 & \sim N(x_1-x_2, 2 \cdot \mu \cdot \Delta t + 2 \cdot \mu \cdot \Delta t) \\ & \sim N(x_1-x_2, 4 \cdot \mu \cdot \Delta t) \;\;\; \text{(21a)} \end{split} $$ $$ \begin{split} Y_1-Y_2 & \sim N(y_1-y_2, 2 \cdot \mu \cdot \Delta t + 2 \cdot \mu \cdot \Delta t) \\ & \sim N(y_1-y_2, 4 \cdot \mu \cdot \Delta t) \;\;\; \text{(21b)} \end{split} $$

Secondly, because the means of the $X_1-X_2$ and $Y_1-Y_2$ are not equal to zero, we need to account for this by calculating the noncentrality parameter $\lambda$. This is equal to the sum of squares of the means of the normalized $X_1-X_2$ and $Y_1-Y_2$ variables, which of course equals (22), corresponding to a scaled version of the original intercellular distance. We then define the standardized variable R’² as (eq. 21). To recap, the meaning of R’² is the scaled squared intercellular distance after updating once with timestep \triangle t starting from an initial distance of R0. Then we define F as the cumulative distribution function (eq. 22). To then finally calculate appropriate timestep for a given r0, rcutoff and µA we take the inverse cumulative distribution function and equate it to the scaled squared intercellular cutoff distance (eq. 23). Note that both the left hand side and right hand side depend on the \triangle t in a nonlinear way, therefore a nonlinear solver is required.

Partial Differential Equations Module

Introductory text about the Partial Differential Equations Module

Computational Scheme

As mentioned earlier the concentrations of AHL and Leucine are modeled using partial differential equations. In the colony level model these equations are solved explicitly. Explicit schemes do not require a lot of work per time step, but unfortunately are not unconditionally stable. In two dimensions the grid ratios $dt/dx^2$ and $dt/dy^2$ can not exceed $dt/dx^2 + dt/dy^2 \leq \frac{1}{2}$ for the solver to be stable. When computing the solution of the hybrid model this requrement forces us to spend a lot of CPU time solving partial differential equations that could be better spent simulation the agents. Therefore an implicit ADI Alternating direction implicit scheme has been implemented. ADI-schemes are unconditionally stable, which allows it to take large time steps with the PDE solver. We used the following scheme: $$ (1 - \frac{1}{2} \mu_x \delta_x^2) U^{n+\frac{1}{2}} + \frac{1}{4}kU^{n+\frac{1}{2}} = (1 + \frac{1}{2} \mu_y \delta_y^2) U^n - \frac{1}{4}kU^n + \frac{\alpha}{2} \rho_A $$ $$ (1 - \frac{1}{2} \mu_y \delta_y^2) U^{n+1} + \frac{1}{4}kU^{n+1} = (1 + \frac{1}{2} \mu_x \delta_x^2)U^{n+\frac{1}{2}} - \frac{1}{4}kU^{n+\frac{1}{2}} + \frac{\alpha}{2} \rho_A $$ In the equations above $\mu$ denotes grid ratios and $\delta^2$ central differences. The production and degradation terms have been incorporated at every time level with a factor of $\frac{1}{4}$. The image below shows the computational molecule of the ADI scheme we chose to implement:

Epanechnikov kernel with h=1

Figure 4
ADI-Molecule. Click to enlarge.

Boundary Conditions

Text about boundary conditions.

Matching

Introductory text about matching

Decoupling of Timesteps

In order to benefit from the implicit PDE-solver described above the agent's time steps are chosen smaller then the time steps of the PDE solver. However type A cells produce molecules continuously as they move trough space. Therefore we record their positions and average over the kernel functions centered at past positions since the last PDE evaluation. That way we avoid blurring the results of the PDE solver too much, when the time step of the agents is reduced. Last but not least we want to point out the relationship between kernel bandwidth and the grid on which the PDE is solved. The larger the kernel bandwidth $h$ is chosen the coarser the PDE grid can be without loosing cells between grid points. If the the diameter of the kernel functions is smaller then the distance between PDE grid points it can happen, that the kernel does not overlap with any PDE grid points. In this case a cells contribution to the molecule concentrations at the next time step could be lost. When the PDE grid is widened to save computing time it is thus necessary to increase the Kernel bandwidth as well. However this increase can lead to a situation where the Kernel function covers significantly more space then the actual diameter of a bacterium. In these cases the Kernel function can be interpreted as a probability function for a cells position. However we avoided too large bandwiths and PDE grids by running our 2-D simulations at the Flemish Supercomputer Center (VSC).

Logo Flemish Supercomputer Center

Figure 6
Logo of the Flemish Supercomputer Center. Click to enlarge.

Bandwidth Tuning

Text about bandwidth tuning.

1-D Hybrid Model



The video box above shows one dimensional simulation results for the hybrid model. A constant speed and random step simulation has been computed. We observe that the bacteria form a traveling wave in both cases, which is essential for pattern formation. These results are also similar to what we get from the continuous model, which confirms our results.

2-D Hybrid Model



The videos above show simulation videos computed at the Flemish supercomputing center, for three different initial conditions similar to the ones we used for the colony level model. The first and second condition start from 9 mixed or 5 colonies of both cell types, arranged in a block or star shape. These first two gradually separate in a manner similar to what we would we also saw in the colony level model. The result for random initial data is fundamentally different. As the agent based approach allows for better implementation of adhesion large cell type A bands form. The AHL and Leucine produced by the type A bacteria causes the B type cells to move away leading to a pattern which we could not produce using PDEs alone, this beautifully illustrates the added value of hybrid modeling.

Incorporation of internal model

Up until now, we have largely ignored the inner life of the bacteria. This inner life consists of transcriptional networks and protein kinetics. Instead we assumed that AHL and leucine production is directly proportional to the density of type A cells. This only works in theory, since bacteria will be affected by their surroundings and the way their dynamics react to it. For example bacteria surrounded by a large concentration of AHL, will have more CheZ and will react more on the presence of Leucine. Also bacteria have different histories and will have different levels of transcription factors and different levels of proteins in their plasma. The proteins are not directly degraded and will still be present in the cytoplasm of the bacteria long after the network has been deactivated. From this, it is clear that 2 bacteria, although surrounded by the same AHL and leucine concentrations, can show different behavior and reaction kinetics.


This results in a heterogeneity of the bacterial population that has not yet been accounted for. To make up for this anomaly, we decided to add an internal model to every agent. This way we will get more realistic simulations. Every agent will get their own levels of CheZ, LuxR, LuxI and so on and will have individual reactions on their surroundings. We hope that this way we can get closer to the behavior of real bacteria.

References

[1] Benjamin Franz and Radek Erban. Hybrid modelling of individual movement and collective behaviour. Lecture Notes in Mathematics, 2071:129-157, 2013. [ .pdf ]
[2] Zaiyi Guo, Peter M A Sloot, and Joc Cing Tay. A hybrid agent-based approach for modeling microbiological systems. Journal of Theoretical Biology, 255(2):163-175, 2008. [ DOI ]
[3] E F Keller and L A Segel. Traveling bands of chemotactic bacteria: a theoretical analysis. Journal of theoretical biology, 30(2):235-248, 1971. [ DOI ]
[4] E. M. Purcell. Life at low Reynolds number, 1977. [ DOI ]
[5] Angela Stevens. The Derivation of Chemotaxis Equations as Limit Dynamics of Moderately Interacting Stochastic Many-Particle Systems, 2000. [ DOI ]
[6] A. L. Schaefer, B. L. Hanzelka, M. R. Parsek, and E. P. Greenberg. Detection, purification, and structural elucidation of the acylhomoserine lactone inducer of Vibrio fischeri luminescence and other related molecules. Bioluminescence and Chemiluminescence, Pt C, 305:288-301, 2000.

Contact

Address: Celestijnenlaan 200G room 00.08 - 3001 Heverlee
Telephone: +32(0)16 32 73 19
Email: igem@chem.kuleuven.be