Team:Valencia UPV/Modeling/diffusion

Since pheromones are chemicals released into the air, we have to consider both the motion of the fluid and the one of the particles suspended in the fluid.

The motion of fluids can be described by the Navier–Stokes equations. But the typical nonlinearity of these equations when there may exist turbulences in the air flow, makes most problems difficult or impossible to solve. Thus, attending to the particles suspended in the fluid, a simpler effective option for pheromone dispersion modeling consists in the assumption of pheromones diffusive-like behavior. That is, pheromones are molecules that can undergo a diffusion process in which the random movement of gas molecules transport the chemical away from its source [1].

There are two ways to introduce the notion of diffusion: either using a phenomenological approach starting with Fick's laws of diffusion and their mathematical consequences, or a physical and atomistic one, by considering the random walk of the diffusing particles [2].

In our case, we decided to model our diffusion process using the Fick's laws. Thus, it is postulated that the flux goes from regions of high concentration to regions of low concentration, with a magnitude that is proportional to the concentration gradient. However, diffusion processes are complex, and modelling them analytically and with accuracy is difficult. Even more when the geometry is not simple (e.g. consider the potential final distribution of our plants in the crop field). For this reason, we decided to consider a simplified model in which pheromone chemicals obey the heat diffusion equation.

Approximation

The diffusion equation is a partial differential equation that describes density dynamics in a material undergoing diffusion. It is also used to describe processes exhibiting diffusive-like behavior, like in our case. The equation is usually written as:

$$\frac{\partial \phi (r,t) }{\partial t} = \nabla • [D(\phi,r) \nabla \phi(r,t)]$$

where $\phi(r, t)$ is the density of the diffusing material at location r and time t, and $D(\phi, r)$ is the collective diffusion coefficient for density $\phi$ at location $r$; and $\nabla$ represents the vector differential operator.

If the diffusion coefficient does not depend on the density then the equation is linear and $D$ is constant. Thus, the equation reduces to the linear differential equation: $$\frac{\partial \phi (r,t) }{\partial t} = D \nabla^2 \phi(r,t)$$

also called the heat equation. Making use of this equation we can write the pheromones chemicals diffusion equation with no wind effect consideration as:

$$\frac{\partial c }{\partial t} = D \nabla^2 C = D \Delta c$$

where c is the pheromone concentration, $\Delta$ is the Laplacian operator, and $D$ is the pheromone diffusion constant in the air.

If we consider the wind, we face a diffusion system with drift, and an advection term is added to the equation above.

$$\frac{\partial c }{\partial t} = D \nabla^2 c - \nabla \cdot (\vec{v} c )$$

where $\vec{v}$ is the average velocity. Thus, $\vec{v}$ would be the velocity of the air flow in or case.

For simplicity, we are not going to consider the third dimension. In $2D$ the equation would be:

$$\frac{\partial c }{\partial t} = D \left(\frac{\partial^2 c }{\partial^2 x} + \frac{\partial^2 c }{\partial^2 y}\right) – \left(v_{x} \cdot \frac{\partial c }{\partial x} + v_{y} \cdot \frac{\partial c }{\partial y} \right) = D \left( c_{xx} + c_{yy}\right) - \left(v_{x} \cdot c_{x} + v_{y} \cdot c_{y}\right) $$

In order to determine a numeric solution for this partial differential equation, the so-called finite difference methods are used. With finite difference methods, partial differential equations are replaced by its approximations as finite differences, resulting in a system of algebraic equations. This is solved at each node $(x_i,y_j,t_k)$. These discrete values describe the temporal and spatial distribution of the particles diffusing.

Although implicit methods are unconditionally stable, so time steps could be larger and make the calculus process faster, the tool we have used to solve our heat equation is the Euler explicit method, for it is the simplest option to approximate spatial derivatives.

The equation gives the new value of the pheromone level in a given node in terms of initial values at that node and its immediate neighbors. Since all these values are known, the process is called explicit.

$$c(t_{k+1}) = c(t_k) + dt \cdot c'(t_k),$$

Now, applying this method for the first case (with no wind consideration) we followed the next steps:

1. Split time $t$ into $n$ slices of equal length dt: $$ \left\{ \begin{array}{c} t_0 &=& 0 \\ t_k &=& k \cdot dt \\ t_n &=& t \end{array} \right. $$

2. Considering the backward difference for the Euler explicit method, the expression that gives the current pheromone level each time step is:

$$c (x, y, t) \approx c (x, y, t - dt ) + dt \cdot c'(x, y, t)$$

3. And now considering the spatial dimension, central differences is applied to the Laplace operator $\Delta$, and backward differences are applied to the vector differential operator $\nabla$ (in 2D and assuming equal steps in x and y directions):

$$c (x, y, t) \approx c (x, y, t - dt ) + dt \left( D \cdot \nabla^2 c (x, y, t) - \nabla \vec{v} c (x, y, t) \right)$$ $$ D \cdot \nabla^2 c (x, y, t) = D \left( c_{xx} + c_{yy}\right) = D \frac{c_{i,j-1} + c_{i,j+1} + c_{i-1,j } + c_{i+1,j} – 4 c_{I,j}}{s} $$ $$ \nabla \vec{v} c (x, y, t) = v_{x} \cdot c_{x} + v_{y} \cdot c_{y} = v_{x} \frac{c_{i,j} – c_{i-1,j}}{h} + v_{y} \frac{c_{i,j} – c_{i,j-1}}{h} $$

With respect to the boundary conditions, they are null since we are considering an open space. Attending to the implementation and simulation of this method, dt must be small enough to avoid instability.

References

The Idea

When one observes moths behavior, they apparently move with erratic flight paths. This is possibly to avoid predators. This random flight is modified by the presence of sex pheromones. Since these are pheromones released by females in order to attract an individual of the opposite sex, it makes sense that males respond to gradients of sex pheromone concentration, being attracted towards the source. As soon as a flying male randomly enters into a conical pheromone-effective sphere of sex pheromone released by a virgin female, the male begins to seek the female following a zigzag way. The male approaches the female, and finally copulates with her [1].

Approximation

In Sexy Plant we approximate the resulting moth movement as a vectorial combination of a gradient vector and a random vector. The magnitude of the gradient vector depends on the change in the pheromone concentration level between points separated by a differential stretch in space. More precisely, the gradient points in the direction of the greatest rate of increase of the function, and its magnitude is the slope of the graph in that direction. The random vector is constrained in this ‘moth response’ model by a fixed angle upper bound, assuming that the turning movement is relatively continuous. For example, one can asume that the moth cannot turn 180 degrees from one time instant to the next.

Our synthetic plants are supposed to release enough sexual pheromone so as to be able to saturate moth perception. In this sense the resulting moth vector movement will depend ultimately on the pheromone concentration levels in the field and the moth ability to follow better or worse the gradient of sex pheromone concentration.

The three clases of male moth behavior we consider for the characterization of males moth behavior are described in Table 1.

This ensemble of behaviors can be translated into a sum of vectors in which the random vector has constant module and changing direction within a range, whereas the module of the gradient vector is a function of the gradient in the field. The question now is how do we include the saturation effect in the resulting moth shift vector. With this in mind, and focusing on the implementation process, our approach consists on the following:

To model chemoattraction, the gradient vector will be always have fixed unit magnitude, and its direction is that of the greatest rate of increase of the pheromone concentration.

To model the random flight, instead of using a random direction vector with constant module, we consider a random turning angle starting from the gradient vector direction.

Thus, how do we include the saturation effect in the resulting moth shift vector? This is key to achieve sexual confusion. Our answer: the behaviour dependence on the moth saturation level --in turn related to the pheromone concentration in the field-- will be included in the random turning angle.

This random turning angle will not follow a uniform distribution, but a Poisson distribution in which the mean is zero (no angle detour from the gradient vector direction) and the standard-deviation will be inversely proportional to the intensity of the gradient of sex pheromone concentration in the field. This approach leads to ‘sexual confusion’ of the insect as the field homogeneity increases. This is because the direction of displacement of the moth will equal the gradient direction with certain probability which depends on how saturated it is.

References

Using a modeling environment called Netlogo, we simulate the approximate moth population behavior when the pheromone diffusion process take place.

The Netlogo simulator can be found in its website at Northwestern University. To download the source file of our Sexy plant simulation in Netlogo click here: sexyplants.nlogo

Setup

• We consider three agents: male and female moths, and sexy plants.
• We have two kinds of sexual pheromone emission sources: female moths and sexyplants.
• Our scenario is an open crop field where sexy plants are intercropped, and moths fly following different patterns depending on its sex.

Females, apart from emitting sexual pheromones, move following erratic random flight paths. After mating, females do not emit pheromones for a period of 2 hours.

Males also move randomly while they are under its detection threshold. But when they detect a certain pheromone concentration, they start to follow the pheromone concentration gradients until its saturation threshold is reached.

Sexy plants act as continuously- emitting sources, and their activity is regulated by a Switch.

The pheromone diffusion process, it is simulated in Netlogo by implementing the Euler explicit method.

Runs

When sexy plants are switched-off, males move randomly until they detect pheromone traces from females. In that case they follow them.

When sexy plants are switched-on, the pheromone starts to diffuse from them, rising up the concentration levels in the field. At first, sexy plants have the effect of acting as pheromone traps on the male moths.

As the concentration rises in the field, it becomes more homogeneous. Remember that the random turning angle of the insect follows a Poisson distribution, in which the standard-deviation is inversely proportional to the intensity of the gradient. Thus, the probability of the insect to take a bigger detour from the faced gradient vector direction is higher. This means that it is less able to follow pheromone concentration gradients, so sexual confusion is induced.

Parameters

The parameters of this model are not as well-characterized as we expected at first. Finding the accurate values of these parameters is not a trivial task. In the literature it is difficult to find a number experimentally obtained. So we decided to take an inverse engineering approach. Doing a model parameters swept, we simulated many possible scenarios, and then we came up with values of parameters corresponding to our desired goal: insects get confused. This will be useful to know the limitations of our system and to help to decide the final distribution of our plants in the crop field. The parameters we found are:

Diffusion coefficient

Range of search: 0.01-0.2 cm^2/s
References: [1], [2], [3], [5]

Release rate (female)

Range of search: 0.02-1 µg/h
References: [4], [5]

Release rate (Sexy Plant)

The range of search that we have considered is a little wider than the one for the release rate of females.
References: It generally has been found that pheromone dispensers releasing the chemicals above a certain emission rate will catch fewer males. The optimum release rate or dispenser load for trap catch varies greatly among species [4]. This emission rate above which male start to get confused could be the release rate from females.

Detection threshold

Range of search: 0.001-1 [Mass]/[Distance]^2

Saturation threshold

Range of search: 1-5[Mass]/[ Distance]^2

Moth sensitivity

This is a parameter referred to the capability of the insect to detect changes in pheromone concentration in the patch it is located and the neighbor patch. When the field becomes more homogeneous, an insect with higher sensitivity will be more able to follow the gradients.
Range: 0-0.0009
(The maximum level of moth sensitivity has to be less than the minimum level of release rate of females, since this parameter is obtained from the difference)

Wind force

Range: -0.1 – 0.1 cm/sec
References: [7] (700cm/sec!! we did not rely on this figure)

Population

The number of males and females can be selected by the observer.

Patches

One can modify the number of patches that conform the field so as to analyze its own case. In our case we used a field of 50x50 patches.

References

The aim consists of reducing the possibility of meeting among moths of opposite sex. Thus, we will analyze the number of meetings in the three following cases:

1. When sexyplants are switched-off and males only interact with females.
2. When sexyplants are switched-on and have an effect of trapping males.
3. When sexyplants are swiched-on and males get confused when the concentration of pheromone level is higher than their saturation threshold.

It is also interesting to analyze a fourth case, what does it happen if females wouldn’t emit pheromones and males just move randomly through the field? :

4. Males and females move randomly. How much would our results differ from the rest of cases?

What is important is that between the first and the third case, the number of meetings should be less in the latter than in the former. Then we are closer to our objective fulfillment.