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

Line 350: Line 350:
 
             <div class="whiterow"></div>
 
             <div class="whiterow"></div>
  
 +
<!--
 
             <div class="center">
 
             <div class="center">
 
               <div class="togglebar">
 
               <div class="togglebar">
Line 747: Line 748:
 
</div><!-- togglebar -->
 
</div><!-- togglebar -->
 
</div><!-- center -->
 
</div><!-- center -->
 +
-->
 +
            <div class="center">
 +
              <div class="togglebar">
 +
                <div class="togglefour">
 +
                    <h2>
 +
                        Hybrid Model
 +
                    </h2>
 +
                </div>
 +
              <div id="togglefour">
 +
                    <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>
 +
            </div>
 +
        <div class="togglebar">
 +
            <div class="togglefive">
 +
                    <h2>
 +
                        Partial Differential Equations
 +
                    </h2>
 +
            </div>
 +
            <div id="togglefive">
 +
                    <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).
 +
                        $$\frac{\partial C(\vec{r},t)}{\partial t}=\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 like most organic
 +
                        molecules 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 (citation).
 +
                        $$ \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}=\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, 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>
 +
              </div>
 +
 +
    <div class="togglebar">
 +
        <div class="togglesix">
 +
 +
                    <h2>Agent-based</h2>
 +
</div>
 +
<div id="togglesix">
 +
                    <p>
 +
                        To model bacteria movement on the other hand, we used an agent-based model that
 +
                        explicitly stored individual bacteria as agents. Spatial coordinates were
 +
                        associated with each agent, which specified their location. After solving the
 +
                        equation of motion for all agents based on their environment, these coordinates
 +
                        were updated at every timestep. In principle, Newton’s second law of motion had
 +
                        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 a constant.
 +
                    </p>
 +
                    <p>
 +
                    <b>
 +
                        Stochastic Differential Equation
 +
                    </b><br/>
 +
                    </p>
 +
                   
 +
 +
 +
                    <p>
 +
                    $$ 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} \;\;\; (6) $$
 +
                        However, the physical “chemotactic force” that propel bacteria is not easily
 +
                        measured or derived. Therefore, we base the equation of motion in one dimension
 +
                        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$,
 +
                        more precisely the (7) PDE for chemotaxis
 +
                        towards some chemoattractant S.
 +
                        $$ \frac{\partial n(\vec{r},t)}{\partial t}=D_n \cdot \nabla^2 n - \nabla (n \cdot \chi(S(\vec{r},t))
 +
                        \cdot\nabla S(\vec{r},t)) \;\;\; (7) $$
 +
That means that when infinitely many particles
 +
                        obey (6), they will exhibit Keller-Segel type spatial dynamics. 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 Ito 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 A$. The first term in (6) on the other hand is easily understood
 +
                        as an advective or drift term (net motion) depending on S, pushing the agents
 +
                        along a positive gradient (for negative chemotaxis the sign is reversed). The
 +
                        chemotactic force 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 defined the chemotactic sensitivity
 +
                        function as in (8).
 +
                        $$ \chi(S(\vec{r},t))=\mu \cdot \frac{\kappa}{S(\vec{r},t)} \;\;\; (8) $$
 +
                        The first important thing to note is that we assumed
 +
                        $\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 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.
 +
                    </p>
 +
                    <p>
 +
                    <b>
 +
                        Cell-cell Interactions
 +
                    </b><br/>
 +
 +
                        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 repulsive-attractive
 +
                        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. 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 bacteria are stochastic and could randomly jump beyond
 +
                        the potential wall, where the force is ill defined. Practically, we’ve decided
 +
                        to define a piecewise quadratic potential (9a),
 +
                        $$ E_{p,attraction}(r_{ij})=\left\{\begin{matrix}
 +
                          0 & 2\cdot r_{cutoff}\leq r_{ij}\\
 +
                          -\frac{1}{2}\cdot k_3 \cdot(r_{ij}-2\cdot r_0)^2 & 2\cdot r_0 \leq r_{ij} < 2 \cdot r_{cutoff} \\
 +
                          -\frac{1}{2}\cdot k_2 \cdot(r_{ij}-2\cdot r_0)^2 &  r_0 \leq r_{ij} < 2\cdot r_0\\
 +
                          -\frac{1}{2}\cdot k_1 \cdot(r_{ij}-\frac{k_1+k_2}{k_1}\cdot r_0)^2 &  0 \leq r_{ij} < r_0
 +
                          \end{matrix}\right. $$
 +
                        $$ E_{p,repulsion}(r_{ij})=\left\{\begin{matrix}
 +
                          0 & 2\cdot r_0\leq r_{ij}\\
 +
                        -\frac{1}{2}\cdot k_2 \cdot(r_{ij}-2\cdot r_0)^2 &  r_0 \leq r_{ij} < 2\cdot r_0\\
 +
                        -\frac{1}{2}\cdot k_1 \cdot(r_{ij}-\frac{k_1+k_2}{k_1}\cdot r_0)^2 &  0 \leq r_{ij} < r_0
 +
                        \end{matrix}\right. \;\;\; (9a) $$
 +
                        $$ \vec{F}_{ij,attraction}(r_{ij})=\left\{\begin{matrix}
 +
                          \vec{0} & 2\cdot r_{cutoff}\leq r_{ij}\\
 +
                          k_3 \cdot(r_{ij}-2\cdot r_0) \cdot \vec{e}_{ij} & 2\cdot r_0 \leq r_{ij} < 2 \cdot r_{cutoff} \\
 +
                          k_2 \cdot(r_{ij}-2\cdot r_0) \cdot \vec{e}_{ij} &  r_0 \leq r_{ij} < 2\cdot r_0\\
 +
                          k_1 \cdot(r_{ij}-\frac{k_1+k_2}{k_1}\cdot r_0) \cdot \vec{e}_{ij} &  0 \leq r_{ij} < r_0
 +
                          \end{matrix}\right.
 +
                        $$
 +
                        $$ \vec{F}_{ij,repulsion}(r_{ij})=\left\{\begin{matrix}
 +
                          \vec{0} & 2\cdot r_0\leq r_{ij}\\
 +
                          k_2 \cdot(r_{ij}-2\cdot r_0) \cdot \vec{e}_{ij}&  r_0 \leq r_{ij} < 2\cdot r_0\\
 +
                          k_1 \cdot(r_{ij}-\frac{k_1+k_2}{k_1}\cdot r_0) \cdot \vec{e}_{ij}&  0 \leq r_{ij} < r_0
 +
                          \end{matrix}\right. \;\;\; (9b)
 +
                        $$
 +
                        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 (9b).
 +
                        Between A type cells, there is a region of attraction (2*r0 < r < 2*rcutoff),
 +
                        where the force points towards the other cell, hence moving them closer together.
 +
                        In the repulsive domain (r < 2*r0), two regions were defined, emulating lower
 +
                        repulsive forces (r0 < r < 2*r0) and higher repulsive forces due to a higher spring
 +
                        constant when the cells are even closer (r < r0). For the purely repulsive interaction
 +
                        scheme there is no attraction and therefore the spring constant for r > 2*r0 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>
 +
                    <p>
 +
                    <b>
 +
                        Equation of Motion
 +
                    </b><br/>
 +
                        Now we are ready to construct the equation of motion for cell type A and B as a
 +
                        superposition of the Keller-Segel SDE (6) and the cell interaction forces,
 +
                        yielding (10).
 +
                        $$  d\vec{r}_{A_i}(t)= \sqrt{2 \cdot \mu_A}\cdot d\vec{W} + \frac{1}{\gamma}\cdot\Bigg( \sum^{A
 +
                            \backslash \{ A_i\}}_j \frac{dE_{p,attraction}(r_{ij}(t))}{dr_{ij}}\cdot \vec{e}_{ij}+\sum^{B}_j
 +
                            \frac{dE_{p,repulsion}(r_{ij}(t))}{dr_{ij}}\cdot \vec{e}_{ij} \Bigg)\cdot dt \;\;\; (10a)
 +
                        $$
 +
                        $$
 +
                            d\vec{r}_{B_i}(t)= \chi(L(\vec{r},t),H(\vec{r},t)) \cdot \nabla L(\vec{r},t)\cdot dt + \sqrt{2
 +
                            \cdot \mu_B(H(\vec{r},t))}\cdot d\vec{W} +
 +
                        $$
 +
                        $$
 +
                            \frac{1}{\gamma}\cdot\Bigg( \sum^{A\cup B\backslash
 +
                            \{ B_i\}}_j \frac{dE_{p,repulsion}(r_{ij}(t))}{dr_{ij}}\cdot \vec{e}_{ij} \Bigg)\cdot dt
 +
                        $$
 +
                        $$
 +
                            \chi(L(\vec{r},t),H(\vec{r},t))= \mu_{B}(H(\vec{r},t)) \cdot \frac{\kappa}{L(\vec{r},t)}
 +
                        $$
 +
                        $$
 +
                            \mu_A(H(\vec{r},t))=\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.
 +
                        $$
 +
                        $$
 +
                            \mu_B(H(\vec{r},t))=\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{(10b)}
 +
                        $$
 +
 +
                        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 has a
 +
                        negative sign signifying that B type cells are repelled by leucine. The 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 switches
 +
                        based on the local concentration of AHL relative to a threshold AHL value, which
 +
                        simulates the dependency of cellular motility on AHL. 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 class="whiterow"></div>
 +
 +
                <div class="widebox">
 +
                  <h2> Heat equation </h2>
 +
                  <p>
 +
                    The left-hand side of (1) is the rate of accumulation of a chemical and the right-hand side is the second spatial derivative of its concentration field. The equation can be understood by considering a one-dimensional concentration profile: 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 to the other side and a result there is no 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 euqation"><img class="example-image" src="https://static.igem.org/mediawiki/2015/0/0a/KU_Leuven_heatDiagram.png" alt="Illustration of the heat euqation" width="45%" height="45%"></a>
 +
<h4><div id=figure1>Figure 1</div> Illustration of the heat euqation. Click to enlarge </h4>
 +
</div><!-- image1 -->
 +
</div><!-- center -->
 +
                    <br/>
 +
                    <!-- 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>
 +
            <!-- video end -->
 +
<div class="whiterow"></div>
 +
</div><!-- widebox -->
 +
               
 +
         
 +
<div class="whiterow"></div>
 +
</div><!-- togglesix -->
 +
</div><!-- togglebar -->
 +
<div class="togglebar">
 +
<div class="toggleseven">
 +
                        <h2>
 +
                            Coupling
 +
                        </h2>
 +
</div>
 +
<div id="toggleseven">
 +
                    <p>
 +
                    <b>Agent-based to PDE</b><br/>
 +
                        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 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).
 +
                        <br/>
 +
                    </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 below. </p>
 +
                          <div class="center">
 +
                        <div id="image1">
 +
                        <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 sum"><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=figure1>Figure 2</div> Kernel sum. 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 an agent in our case.
 +
                        Some of the most common kernel functions include gaussian kernels, triangular
 +
                        kernels and epanechnikov kernels.
 +
                        During our simulations we have found the epanechnikov kernel particularly useful.
 +
                        In two dimensions these are defined as:
 +
                        $$k(x,y) = \left\{\begin{matrix} \frac{3}{4h} * (1 - ((x/h)^2 + (y/h)^2)) & \text{if } ((x/h)^2 +
 +
                                  (y/h)^2) \leq 1 \\
 +
                                  0 & \text{else}  \end{matrix}\right. \;\;\; (11) $$
 +
                        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. 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. Using a KDE allows us to define the agent density at any
 +
                        point by referring to the kernel sum.
 +
                        </p>
 +
 +
                      <div class="center">
 +
                        <div id="image3">
 +
                        <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 with h=1"><img class="example-image" src="https://static.igem.org/mediawiki/2015/5/52/KU_Leuven_EpanechnikovDifferentBandwidth.png" alt="Epanechnikov kernel with h=1" width="45%"></a>
 +
                        <h4><div id=figure3>Figure 3</div> Epanechnikov kernels with bandwidths with various bandwidths. Click to enlarge </h4>
 +
                        <a class="example-image-link" href="https://static.igem.org/mediawiki/2015/8/88/KU_Leuven_variousKernelTypes.png" data-lightbox="example-set" data-title="Epanechnikov kernel with h=2"><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 4</div> Gaussian, triangular and Epanechnikov kernel functions. Click to enlarge </h4>
 +
</div><!-- image3 -->
 +
                <div class="whiterow"></div>
 +
                  <p>               
 +
                  <b> PDE to agent-based </b><br/>
 +
 +
                  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, which means that they can be located
 +
                  at any point of the domain.
 +
                  As we have seen, 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. 
 +
                </p>
 +
</div><!-- center -->
 +
</div><!-- toggleseven -->
 +
</div><!-- togglebar -->
 +
</div><!-- center -->
 +
 
<div class="whiterow"></div>
 
<div class="whiterow"></div>
 
<!--</div>-->
 
<!--</div>-->

Revision as of 14:30, 9 October 2015

Logo

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 modelling 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

Hybrid model

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}=\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 (citation). $$ \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}=\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, 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.

Agents

To model bacteria movement on the other hand, we used an agent-based model that explicitly stored individual bacteria as agents. Spatial coordinates were associated with each agent, which specified their location. After solving the equation of motion for all agents based on their environment, these coordinates were updated at every timestep. In principle, Newton’s second law of motion had 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 a constant.