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

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$.

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).

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.

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 4).

Figure 4

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 5. 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.

Figure 5

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

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.

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).

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.

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 6.

Figure 6

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.

Figure 7

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

In practice, scaled versions of standard kernel functions are used, which are of the form (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.

Figure 8

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.

Figure 9

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.

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.

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.

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 (11) 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 executed 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 one million. This puts a heavy computational load on the model which can be reduced greatly by using better algorithms. The crucial observation here is that the force is zero when the distance between two cells is larger than $2\cdot r_{cutoff}$. An agent therefore only actually interacts with its nearest neighbors. Hence, the naive implementation wastes a lot of time computing interactions which have no effect.

To speed up calculation of cell interactions we used an algorithm that searches only for neighbors within the cutoff distance for every agent. Moreover, since we take the timestep small enough so that agents do not easily move out of each other's influence from one moment to another (see paragraph below), it is not necessary to recompute the neighbors at every timestep. Instead, we maintain a neighbor list for every agent which is refreshed only every few timesteps.

There exist many different algorithms for searching nearest neighbors and the optimal choice of algorithm depends on the characteristics of the spatial distribution of agents. Since in our model cells 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. Determining which cell an agent belongs to is easily done by rounding the x and y-coordinates down to the nearest multiple of the grid spacing. In Figure 10 agent 1 and agent 2 have different coordinates but $x \; mod \; d = x_k$ and $y \; mod \; d = y_k$ for both agents and therefore they will be assigned to the same cell k.

Figure 13

Cell assignment. Click to enlarge.

If the grid spacing is chosen so that interactions do not reach further than the adjacent cells, every agent only needs to look for neighbors in 9 cells (its own cell plus 8 adjacent cells) instead of the entire domain, as illustrated in Figure 12.

Figure 14

Neighbor search in adjacent cells. Click to enlarge.

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).

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).

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$, 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).

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).

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).

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 $\bar{R}^2$ as (23).

To recap, the meaning of $\bar{R}^2$ is the scaled squared intercellular distance after updating once with timestep $\Delta t$ starting from an initial distance of $r_0$. Then we define $F$ as the cumulative distribution function for the random variable $\bar{R}^2$ as (24).

To then finally calculate appropriate timestep for a given $r_0$, $r_{cutoff}$, $\mu$ and desired $P^-$ we take the inverse cumulative distribution function and equate it to the intercellular cutoff distance, after appropriate scaling and squaring (25).

Note that both the left hand side and right hand side both depend on the $\Delta t$ in a nonlinear way, therefore a nonlinear solver is required to solve this equation. As stated previously, it is not possible to completely avoid that the cells will jump out of each other's influence in a single timestep. However, equation (25) allows us to calculate the timestep $\Delta t$ that will result in cells detaching in a single timestep with a probability $P^-$ that we can specify.

Solving the equation of motion to compute trajectories of bacteria is well-defined as long as they stay within in the domain. To give a complete picture however, we also need to specify what agents should do when they attempt to cross the spatial boundaries of the simulation. Initially we implemented a simple “freeze” boundary condition, which means that the bacteria are simply forced back to the edge whenever they tried to exit the domain. However, as our model got more sophisticated it became apparent that simulating an entire petri dish would not be feasible. Therefore, under the assumption that our system would produce similar patterns over the entire petri dish, we decided to implement periodic boundary conditions. This means that we treat the simulated domain as a “unit cell” which is repeated ad infinitum in all spatial directions. As a result, an agent that tries to exit the simulation box at one side re-enters from the opposite side, like in Figure 13.

Figure 15

Periodic boundary conditions. Click to enlarge.

Periodic boundary conditions are easily implemented by taking the modulus of every coordinate with the length of the periodic domain in that dimension as modulo (26).

One last complication is that cell-cell interactions across domain edges are possible with periodic boundary conditions. To compute those, virtual cells containing virtual agents are used that represent real cells and agents close to the opposite edge of the domain, but are translated to appear to come from outside the domain.

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 timestep, but unfortunately are not unconditionally stable. If $\Delta t$ the timestep, $\Delta x$ and $\Delta y$ the distances between grid points in the x and y-direction respectively, and $D$ the diffusion coefficient, the condition for stability is given by (27).

When computing the solution of the hybrid model this requirement forces us to spend a lot of CPU time solving partial differential equations that could be better spent simulation the agents. Therefore an implicit alternating direction implicit (ADI) scheme was implemented. ADI schemes are unconditionally stable, which allows us to take large time steps with the PDE solver. The difference from a To derive the scheme we first reiterate the partial differential equation (4).

We could use a fully implicit scheme by taking the field values on the right-hand side at $t+\Delta t$, but there is no way to reshape the resulting matrix to take on a banded form and hence no efficient solver is applicable to compute the field values at the next timestep. The trick used in ADI is then to divide the timestep $\Delta t$ in two half timesteps and take only one dimension implicitly at a time, and leave the other dimension in its explicit form. In the subsequent half timestep, the dimensions are switched (hence alternating). For example, if we take $x$ implicitly in the first half, we obtain (28a).

The density of type A bacteria $\rho_A$ might evolve over multiple timesteps in $[t;t+\Delta t]$ due to decoupling of timesteps (see Matching below), therefore we compute an average over the entire interval and use that in the discretization scheme. Furthermore, the degradation term can be incorporated in a Crank-Nicholson manner since they do not depend on neighboring grid points and therefore do not influence the shape of the matrix. We define $r_x$ and $r_y$ as in (29).

Inserting these in (28a) and moving all implicit terms to the left-hand side and all explicit terms to the right-hand side gives (28b).

The next half timestep is analogous, but with the finite difference in the y-direction taken implicitly.

In the end we obtain a unconditionally stable scheme to calculate the evolution of AHL and leucine concentrations over the domain. Moreover, this scheme can be computed efficiently by solving banded matrix systems. The image below illustrates the computational molecule of the ADI scheme we implemented:

Figure 16

ADI-Molecule. Click to enlarge.

As explained above in the implementation of the agent-based module, we chose to use periodic boundary conditions for our system. The field should thus repeat itself in all directions. In the PDE-based module, this comes down to equating the “missing” neighbor of a grid point on the edge of the domain to the grid point at the opposite edge of the domain. Assuming our grid consists of $J_x\;x\;J_y$ points, the grid point to the “left” of $C_1,j$ is for example $C_{J_x},j$. Conversely the grid point to the “right” of $C_{J_x},j$ is $C_1,j$. This impacts our computational scheme through a modified finite difference expression for grid points on the edge of the domain. The regular form is defined as (31).

However, for the outermost grid points of our domain the finite difference is defined as (32).

Implementing these adjusted finite difference expressions in our computational schemes breaks the bandedness of our matrices, but the result is still a nearly-banded sparse matrix for which efficient methods also exist. Therefore, we can implement periodic boundary conditions efficiently.

In order to benefit from the implicit PDE-solver described above the agent's timestep is chosen smaller then the timestep of the PDE solver. As a result, the agent density field goes through several timesteps before the PDE solver takes one timestep. The question then becomes at which timestep the density field should be sampled to provide the source term in the PDEs. If only the density field at one particular timestep is chosen to represent the density field during the entire interval $[t;t+\Delta t]$, it is clear that a lot of information is lost because type A cells produce molecules continuously as they move trough space.

Therefore we record the density fields at all timesteps and since the last PDE evaluation and average of them. That way we avoid skewing the results of the PDE solver too much when the time step of the agents is reduced. In the videoplayer below we illustrate this using a simple system involving a single agent moving along a quasi-1-dimensional domain with a smaller timestep than the PDE solver. When the density is not averaged, the production of AHL and leucine appear to be spiked at regular intervals. However, when our averaging method is used the production shows a more uniform profile, which is closer to the ideal case of continuous production.

As we scaled up our model to realistic dimensions we found that the discretization of the PDE grid had to grow very large relative to the size of bacteria. Therefore, if we kept the bandwidth constant while scaling up the grid size, it would become possible for the grid to "miss" the bacteria. Indeed, as the reach of the kernel functions stays constant and the distance between adjacent grid points becomes larger than the diameter of kernel functions, agents will inevitably start falling off the radar. In the extreme case of small bandwidth and large grid spacing, no agents will be picked up by density calculations at all and their contribution to AHL and leucine production will be completely lost. Therefore, as the PDE grid is scaled up, the bandwidth needs to be increased accordingly. However this inflation of the bandwidth can lead to a situation where the kernel function covers significantly more area than what the actual bacterium physically occupies. In this case the interpretation the kernel function can be altered to represent a probability function for a cell's position. We were able to avoid increasing the bandwidth to unrealistic proportions while scaling up the PDE grid dimensions by running our simulations at the Flemish Supercomputer Center (VSC).

Figure 17

Logo of the Flemish Supercomputer Center. Click to enlarge.