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 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 of the
simulations are shown and summarized. Finally, the incorporation
of the internal model into the hybrid model is discussed.
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.
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.
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.
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.
Implementation
Agent-Based Module
In this paragraph, we will discuss the implementation of the agent-based module. First, the problem of computing pairwise interactions is discussed and a nearest-neighbor algorithm is introduced to reduce the computation time. Then, the choice of timestep is considered with regards to the cell-cell interaction scheme. Finally, the boundary conditions for agents are defined and explained.
Partial Differential Equations Module
To implement the partial differential equations (4) we have to discretize our domain into a grid and apply some computational scheme that describes the evolution of our fields in discrete timesteps. We will discuss only the scheme that we employed for the two-dimensional model. Furthermore, boundary conditions are described below as well.
Matching
The purpose of the hybrid model is to employ appropriate modeling techniques for different entities that behave differently and couple them within a single system. The advantage is that computational costs can be lowered while retaining accuracy. However, if no special care is taken to adapt the modules to each other in an efficient way, the computational savings might not be significant or even nonexistent. Therefore in this paragraph we discuss the matching of temporal scales and spatial scales in the hybrid model.
Results
The videos above show the results of simulations that were run
at the Flemish Supercomputer Center.
Three different initial conditions similar to those used in
the colony level model were used.
The first and second initial conditions correspond to
colonies of both cell types arranged into a block and
a star-shaped pattern respectively. Cell type B can be seen
fleeing the colonies due to the production of AHL and leucine
by cell type A. Therefore the system evolves qualitatively in
the same way as in the colony level model.
However, the model exhibits very different behavior when
the initial positions of bacteria are random. Whereas
the colony level model did not evolve into any convincing
pattern, cell types A and B clearly separate into different
domains in our hybrid model. Cells of type A cluster into blob-like
formations while cells of type B fill up the space in between.
Our hypothesis therefore is that the sticky membrane proteins
expressed by cells A is crucial for the formation of regular patterns.
To verify this, we performed a negative control by disabling
the attractive cell-cell interactions and running the simulation
under the same conditions. The results are similar to the
colony level model in the sense that there is very little to no pattern
formation to be seen. Therefore we can conclude from the hybrid model
simulations
that cell-cell interactions will play a crucial role in the pattern
formation in our system.
Furthermore, by developing the hybrid model in parallel with
the colony level model we were able to compare the results and
demonstrate the
usefulness of incorporating an agent-based module into our
modeling approach. It allowed us to calculate cell-cell interactions
into our simulations which were otherwise not possible to model
and eventually lead us to valuable
insights into the mechanisms of pattern formation.
Incorporation of Internal Model
In this model, we have ignored the inner life of bacteria. This includes gene expression, signal transduction, quorum sensing etcetera. Instead we assumed that AHL and leucine production is directly proportional to the density of type A cells. As the internal model however demonstrates, this is a very rough approximation. Bacteria behave differently based on their surroundings and the signals they have been exposed to earlier. After all, protein concentrations can fluctuate due to external disturbances, so the response to signals from the environment can be heterogeneous across a bacterial population. We can simulate this dynamic response using the internal model and therefore incorporating the internal model into the hybrid model provides a path for further improvement of our model and possibly more clues into the mechanisms of pattern formation in nature. However, during our project we have decided to focus our resources on the individual models therefore we leave this next step for future iGEM teams to explore.
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. |