We decided to model the deployment of our bed bug trap in a realistic environment, to understand how it was likely to work and what parameters might be important. This task involved a number of different modeling tasks, that posed different problems: modeling the diffusion of pheromones in the room coming both from a natural bed bug nest and the traps that we plan to deploy; modeling the movement of bed bugs influenced both by the pheromone field and their nest; and modeling the fungal epidemic that we plan to induce in the bed bug population.

Though the model is complex and includes many parameters, it has already allowed us to draw several conclusions, and as the model is improved and the parameters are refined other conclusions will follow. For instance, using realistic diffusion parameters it is clear that the pheromone field is relatively rapidly established and so the precise nature and concentration of pheromones is probably not critical. In contrast, the delay between infection and death of bed bugs is critical for ensuring the eradication of the nest. Our modeling thus helps understand critical aspects of the proposed design.

How did we model?

Our program has to model the diffusion of pheromones through a room containing furniture. It needed to include the behaviour of bed bugs influenced by the pheromones and to describe their deaths (caused by the fungus). This model will have to manage a lot of parameters: from the diffusion of pheromones to the bed bugs movements and the interactions between all these variables. Because of the complexity of the phenomena, build the program by describing every action occurring in our model would be too difficult. Fortunately for us, we could use the language NetLogo, which allowed us to use agent-based modelling, a type of modelling particularly efficient for complex system modelling.

Agent-based modelling allows us to define the behaviour of the agents (the diffusing molecules, the bedbugs,...) by defining the interactions between those agents. With the NetLogo language, those agents are patches, representing the environment where other agents are interacting. The turtles are mobile agents who can move through the world made by patches and interact with other agents (patches or turtles). Those agents allowed us to build our current program. Thanks to this, it's easier to model systems with multiples agents interacted with every other agents is simpler.

What did we Model ?

Modelling Diffusion

When we begin the development of our project, we first searched on the iGEM's database, in order to find if any teams had already worked on pheromone diffusion. After a few research among the previous iGEM's teams, we found out that the Valencia team of 2014 had also worked on pheromone diffusion and had modeled it using NetLogo. [1]

Because of the similarity between their work and ours, we inspired ourselves with the design of their modelling of diffusion when we model this aspect.

Fick's law of diffusion, mathematical model

As the 2014 Valencia team already found out, there is many ways to model diffusion. Here, we will use the Fick's law of diffusion to model it because of its simplicity and the amount of documentation about it we have access to. But there could be an impact of diffusion within a fluid in motion, the effect of drafts on diffusion. This can be modeled by the Navier-Stokes equations but since the nonlinearity of this equation make it a very complex (if not impossible) equation to solve, we choose to not use it and hypothesize a room without drafts.

How the equation is used in our model

The most well known form of the Fick's law of diffusion is described as it follows :

\begin{equation} \frac{\partial c}{\partial t} = {\nabla(D,c)} \end{equation}

It is a partial differential equation with D being the diffusion coefficient, c being the pheromone concentration and \(\nabla\) being the gradient operator. Since we consider our diffusion coefficient as constant, the equation can be set as:

\begin{equation} \frac{\partial c}{\partial t} = D\times\nabla^2 c = D\times\Delta c \end{equation}

With \(\Delta\) being the Laplacian operator. This equation is called the Heat equation and is well known to model various transfer problems. In our case, this equation models the diffusion of pheromones, in the absence air flow. In a 2D environment, the equation would be :

\begin{equation} \frac{\partial c}{\partial t} = D\times\Delta c = D\times(\frac{\partial^2 c}{\partial x^2} + \frac{\partial^2 c}{\partial y^2}) \end{equation}

Adapt the equation to our model, discretization

Even though we have our mathematical model, we still have to adapt it to our program. In order to do it, we use the finite difference method to approximate the solution of our equation. This method allows the program to compute on every position the approximate solution of our equation. By this, we are able to monitor the diffusion of the phermones at each point in time and space. Finite difference method can be achieved through various ways. Here, and even though the implicit method is stabler we choose the Euler's explicit method, since it is the simplest method for us to use.

This method allow discretization through 3 steps :

• Splitting time and space

We firstly need to split each time an space variable in order to compute the variable at each point and time.

Here, we splitted the coordinates with the same value, noted \(\Delta a\) (with \(a\) either x, y or t) for each type of coordinates. It is written as follows :

\begin{equation} x_i = i\times\Delta x \end{equation} \begin{equation} y_j = j\times\Delta y \end{equation} \begin{equation} t_n = n\times\Delta t \end{equation}

• Discretization of the equation

We now need to translate the equation from a continuous to a discrete way, in order to make the computation of the solution possible. There is 3 elements to be descretized in our equation: \(\frac{\partial c}{\partial t}\), \(\frac{\partial^2 c}{\partial x^2}\) and \(\frac{\partial^2 c}{\partial y^2}\). In the first case, the partial derivative can be discretized as follows:

\begin{equation} \frac{\partial c}{\partial t} = \frac{c^{n+\Delta t}_{i,j} - c^n_{i,j}}{\Delta t} \end{equation}

With \(c(x_i,y_j,t_n) = c^n_{i,j}\), and so \(c(x_i,y_j,t_n + \Delta t) = c^{n + \Delta t}_{i,j}\) For the other 2 partial derivatives, they wil be discretized as follows :

\begin{equation} \frac{\partial^2 c}{\partial x^2} = \frac{c^n_{i + \Delta x,j} + c^n_{ i - \Delta x,j} - 2\times c^n_{i,j}}{\Delta x^2} \end{equation} \begin{equation} \frac{\partial^2 c}{\partial y^2} = \frac{c^n_{i,j + \Delta y} + c^n_{i,j - \Delta y} - 2\times c^n_{i,j}}{\Delta y^2} \end{equation}

The equation can than be written as :

\begin{equation} \frac{c^{n+\Delta t}_{i,j} - c^n_{i,j}}{\Delta t} = D\times(\frac{c^n_{i + \Delta x,j} + c^n_{ i - \Delta x,j} - 2\times c^n_{i,j}}{\Delta x^2} + \frac{c^n_{i,j + \Delta y} + c^n_{i,j - \Delta y} - 2\times c^n_{i,j}}{\Delta y^2}) \end{equation}

• Solving the equation thanks to discretized equation

Once simplified, the approximate solution of the discretized equation can be written as follows :

\begin{equation} c^{n + \Delta t}_{i,j} = D\times(\frac{c^n_{i + \Delta x,j} + c^n_{ i - \Delta x,j} - 2\times c^n_{i,j}}{\Delta x^2} + \frac{c^n_{i,j + \Delta y} + c^n_{i,j - \Delta y} - 2\times c^n_{i,j}}{\Delta y^2})\times\Delta t + c^n_{i,j} \end{equation}

To avoid instability, we need to fulfill the condition \(\frac{D\times\Delta t}{\Delta x^2}\leqslant \frac{1}{2}\).

Modelling the impact of furnitures

Since pheromones will diffuse in a room containing furniture, we had to take into account that they could impact on diffusion. We will integrate it to our program thanks to the computation of the effective diffusion coefficient \(D_eff\) [2]. The computation will be realized thanks to a modified Maxwell-Garnett equation.

This equation is generally used to solve problems of transport in composite materials. Here, we will use it to compute the effective diffusion coefficient in a case where the environment is composed of square inclusions inside the general environment (here the furnitures included inside the room).

Usually the Maxwell-Garnett equation (if the diffusion is studied) is described as follow:

\begin{equation} \frac{D_eff - D_m}{D_eff + 2\times D_m} = \phi(\frac{D_i - D_m}{D_i + 2\times D_m}) \end{equation}

With \(D_m\) being the diffusion coefficient of the environment, \(D_i\) being the diffusion coefficient inside inclusions and \(\phi = \frac{V^d_i}{V^d_m}\) with \(V_i\) the volume of inclusions and \(V_m\) the volume of the environment. Because of our 2D environment, we will compute areas instead of volumes. They will be compute with the length of the furniture (for the area of the inclusions) and the length of the entire world (for the area of the environment). From this equation, and with an environment in 2D, we can set \(D_eff\) as follows:

\begin{equation} D_eff = D_m \times (1 + \frac{2 \times \phi \times (D_i - D_m)}{D_i + D_m - \phi \times (D_i - D_m)}) \end{equation}

But, here, inclusions (furniture) can't be crossed by pheromones. After applying a correction to the equation in order to take account of this impenetrability, the equation describing \(D_eff\) in our case can be written as follow:

\begin{equation} D_eff = \frac{D_m}{1 + \phi} \times (1 - \frac{2 \times \phi}{1 + \phi}) \end{equation}

Modelling the bed bugs behaviour, movements and death

Movements of bed bugs

Now we have modelled the diffusion of pheromones through the room, now, we have to model the behaviour of bed bugs in this environment. We firstly visualised the movement of our bed bugs as random when they are not in contact with any pheromones. Then, as soon as they are in contact with pheromones, the randomness of their behaviour is limited by the attraction of the high quantity of pheromones. Since pheromones are used to attract bed bugs to the traps, only 2 structures will emit pheromones : traps, in order to model our solution, and the nest, modelling the place where the bed bug go back.

Concretely, when the bed bug is not in contact with a patch containing pheromones, we designed the bedbug to walk randomly across the environment, thanks to a group of ifelses instructions and an increasing randomness code. When the bedbug is near a patch containing pheromones, the bed bug is oriented facing the patch with the higher quantity of pheromones. Then, the bed bug moves forward to the patch the bedbug is facing, but with a bit of randomness again. It is important to note that bed bugs have a minimum and a maximum pheromone detection threshold. Those thresholds limit the attractive effect of the pheromones, making the model more realistic.

We also implemented some particular movements behaviour in the case a bedbug had to interact with a specific patch, such as furniture, traps or the nest. When the bedbug is facing a furniture, it will be oriented at the opposite from furniture before moving forward randomly. In the case of the nest, the bed bug is stopped by the histamin contained in the nest. A few moments later, the bed bug is set facing another direction and move forward.

Death mechanisms of bed bugs and their nest

Another particular behaviour is triggered by the infection of bed bugs by the fungus and the deaths of infected bed bugs. When a bed bug is in front of a trap, the program colors the bed bug in orange and set off a timer. It represents the life time left to the bed bug from its infection and this is unique for each bed bug. When that timer is finished, the bed bug dies and colors the patch below in yellow. That patch, the corpse of the bed bug, can also infect other bedbug and contribute to the propagation of the infection.

Another particular situation is when an infected bed bug dies next to the nest. In that case, the fungus will spread to bed bugs inside the nest, adults and younger individuals (like eggs). In order to model this, we give the nest a value, representing the number of time the nest can withstand the death of an infected bed bug next to it. Each time a bed bug dies this way, this counter is decreased by one, and when it gets to 0, the nest is destroyed, representing the infection of every individuals inside the nest.

Simulation routine

The Netlogo environment can be downloaded on it website, at the Nothwestern University [3]. The source file of our program is available here.

Initialization

When the program is set up, 4 types of agents are created:

• bedbugs: each of them will have 2 major features. The lifetime available after the infection, the amount of time the bedbug will be immobilized in front of the nest. They are all created inside the nest.
• traps: they will diffuse pheromones and start the infection of bedbug.
• furnitures: impenetrable obstacles for any other agents
• nest: location where the bedbugs go back after exploration

During this step, the program also compute the effective diffusion coefficient, and prepare the environment for diffusion.

Running Cycle

During a run, the program computes the quantity of pheromones in every patch and colors it accordingly. The more pheromones are in the patch, the redder the patch is. Thanks to this computation, the bedbugs can guide themselves towards the higher quantity of pheromone.

After this calculation step, the program manages the specific interaction between bedbugs and the nest. Following this step, the bedbugs moves through the environment. Finally, the infection and deaths mechanisms are setted up.

Results