Team:KU Leuven/Modeling/Hybrid
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.
