Modeling: Spatial Approach

In this chapter we will present the model we used to do spatial simulations of our interval switch network. This network is intended as part of a cellular automaton that serve as a detector that can differentiate between having ’too few’ and ’too many’ neighbors, in this chapter we will try to demonstrate this concept for the case where no neighbors is ’too few’ and two neighbors is ’too many’.

The network we simulate is shown in Figure 6.7. It consists of one hybrid promoter, inhibited by CI and activated by HSL.LuxR, which facilitates GFP expression. The second promoter is a modified lux promoter which is activated by HSL.LuxR and facilitates CI expression. Because the lux promoter is modified in such a way that the binding site is less easily bound to, we expect that more HSL.LuxR is needed to express enough CI to block the hybrid promoter. This means that when HSL. LuxR is too low, the hybrid promoter (and GFP expression) will not be turned on, and when HSL.LuxR exceeds a certain treshold, the lux promoter will be activated and its CI expression will block the hybrid promoter, stopping GFP expression. The result is a network that will only show a GFP response if the HSL.LuxR concentration is inbetween two limits.

Figure 6.7. Device diagram.

In this section we will give a specification of our model on a ’high level’, meaning that each sensing and protein producing part is modelled as a single Ordinary Differential Equation (ODE) using Hill equations [23]

. For example a Lux-promoter that is connected to a Ribosomal Binding Site (RBS) and is followed by a protein coding region will be modelled using a single Hill equation that returns the protein production as a function of the activators/inhibitors of the promoter. In this way, each of the producing parts that depend on promoter activation or inhibition can be modeled using a single Hill Equation of the form equation 1 for promoters that are activated by Sa and equation 2 for promoters that are inhibited by Si. Where Vmax represents the maximal expression (without inhibition and/ or maximal activation) of the protein and the parameter Km represents the concentration of the inhibitor/ activator at which the production rate is half-maximal.

Of course we also have to model ’standard’ linear reactions such as decay (reaction 48). We do this using Mass Action (MA) kinetics (resulting in equation 49). Note that species depicted in brackets [], represent the concentration.

In the following subsections we will present the model we used to simulate the behavior of our gene network. The Matlab source code itself can be downloaded from 'Modeling Files'

Species

This section lists the different species that are present in the model.

Parameters

The parameters needed for our model have been taken mostly from a single paper [24]

. and estimated from the values found there. Other sources for parameters are [9]

S. Basu, R. Mehreja, S. Thiberge, M. T. Chen, and R. Weiss. Spatiotemporal control of gene expression with pulse-generating networks. PNAS USa, 101(17):6355–6360, April 2004.

and [25]

D. Braun, S. Basu, and R. Weiss. Parameter estimation for two synthetic gene networks: a case study. Acoustics, Speech, and Signal Processing, Proceedings (ICASSP ’05), 5. IEEE International Conference, March 2005

.

In this section we list the Differential Equations that make up our model. Species depicted in brackets [], depict the concentration of a substance. Each differential equation represents the derivative of the concentration of a substance, in other words the production of a species at a single time-step.

These ODEs take into account the diffusion of HSL to and from other cells (the current cell is (x,y)). For this model we consider a 2-dimensional grid and it is assumed that each cell has four neighboring cells.

LuxR

Because LuxR is constitutively expressed, we simply keep the concentration of LuxR constant at 0.5μM. This Same approach is used in [24]

for the wild-type LuxR.

HSL and LuxR have to form a complex first before they can form the HSL.LuxR complex, this is modeled in de ODE by quadratizing both terms [24]

.

GFP response depends on HSL.LuxR and CI concentrations, where the first acts as an activator and the second serves as an inhibitor.

The amount of CI expressed depends on the modified lux promoter, meaning that expression is activated by the amount of HSL.LuxR. Because this promoter is weaker then the promoter that activates GFP expression, i.e. it needs more HSL.LuxR, we can create a band of HSL.LuxR concentrations where GFP is expressed but CI is not yet (sufficiently) expressed.

To solve this collection of ODEs we used MATLAB’s [11]

ode15s function. The .m file we used for the simulations can be found on the website and contains all the listed ODE functions. To study the response of our model we included sender populations, referred to below as sender cells, in our spatial grid. These populations constantly produce 1μM of HSL and are depicted in the images as red cells. Because we are interested in the steady states, the images shown are snapshots taken after running the model for approximately 1000 time-steps (minutes). In the following sections we will explore several parameters and try to approach our desired result, i.e. an interval switch that is tweaked for use in a cellular automaton. The grid-size used in the following experiments is 17 x 17 (’uninteresting regions’ (i.e. without GFP expression) are cropped from the images).

Standard Parameters

When we simply run our model with the parameters described in the paper and include a single sender cell, we get the result as shown in Figure 6.8. We can see that, as expected, a band of GFP forms around the sender cells, indicating that the cells that are near a sender get too much HSL and the cells that are far away from the sender does not get enough HSL for proper activation.

Figure 6.8. GFP response with a single sender cell.

Decreasing lux Promoter Sensitivity

After setting the concentration of HSL.LuxR that is needed for half-maximal expression of CI twice as high (to 0.16μM), we get the result shown in Figure 6.9. There we can see that more HSL is required to produce enough CI to inhibit GFP output, slightly increasing the width of the band.

Figure 6.9. GFP response with a single sender cell, with Kl(a)increased to 0.16μM.

Standard Parameters

When we run the model with the normal parameters and two sender cells, we get the result as shown in Figure 6.10. This indicates that cells that receive HSL from both senders now approach the high threshold and show a lesser GFP response.

Figure 6.10. GFP response with two sender cells (default parameters)

Tweaked Parameters

We now try to adjust the parameters so that we get a response that is suitable for use in a cellular automaton. In this case we try to tweak the model so that when two neighboring cells are ON, a cell does not show a GFP response, and when only a single neighboring cell is ON, a cell does show a GFP response. After playing around with the parameters we found that one way of achieving this behavior is by drastically decreasing HSL diffusion to 0.00001 mm2 min-1 and increasing HSL decay to 0.03 min-1, the resulting snapshot of the grid is shown in Figure 6.11.

Figure 6.11. GFP response with two sender cells. HSL diffusion set to 0.00001mm2 min-1 and HSL decay set to 0.03 min-1.

A nice thing to note is that such parameter-tweaks can also be done in ’real biology’ like by modifying the pH of the solution to increase HSL decay or increasing the density of the medium or distance between cells to reduce HSL diffusion speed.