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 frequent nonlinearity of these equations makes most problems difficult or impossible to solve, since it may exists turbulences in the air flow [1]. Now attending to the particles suspended in the fluid, an 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 [2].

There are two ways to introduce the notion of diffusion: either 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.

In our case, we decided to hold our diffusion process by the Fick's laws. So 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 (final distribution of our plants in the crop field). For this reason, we decided to consider a simplified model in which pheromone chemicals obey to the heat diffusion equation.

Approximation

The diffusion equation is a partial differential equation which 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 following 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 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 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 that the quantity is moving. Thus, $\vec{v}$ would be the velocity of the air flow.

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

For determining a numeric solution for this partial differential equation are used the so-called finite difference methods. The technic consists in approximating differential ratios as $h$ is closer to zero, so they are useful to approximate differential equations.

With finite difference methods, partial differential equations are replaced by its approximations in finite differences, resulting in an algebraic equations system. The algebraic equations system is solved in each node $(x_i,y_j,t_k)$. These discrete values describe the temporal and spatial distribution of the unknown function.

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.

Euler explicit method is the simplest option to approximate spatial derivatives, in which all values are assumed at the beginning of Time.

The equation gives the new value of the pheromone level in terms of initial values in 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 implies that the expression that refers to 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, it is applied central differences to the Laplace operator $\Delta$, and the backward differences 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 opened-space. Attending to the implementation and simulation of this method, dt must be small enough to avoid instability.

The Idea

When one stares at moths, they apparently move with erratic flight paths. It is possibly because of predator avoiding reasons.

In this frame, sex pheromones influence in moth behavior is also considered. 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 comes into conical pheromone-effective sphere of sex pheromone released by a virgin female, the male begins to seek the females in a zigzag way, approaches to the females and finally copulates with her [1].

Approximation

In this project we approximate the resulting moth movement as a vectorial combination of a gradient vector and a random vector. The magnitude of the gradient vector comes from the change in pheromone concentration level among 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 restricted in this ‘moth response’ model by a fixed angle, assuming that the turning movement is relatively continuous and for example the moth can’t turn 180 relative degrees at the next instant.

Since the objective of this project consists in avoiding pest damage by reaching the mating disrupting among moths, our synthetic plants are supposed to release enough sexual pheromone so as 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.

At this point, let’s highlight the three main aspects we consider for the characterization of males moth behavior:

In this context, this ensemble of behaviors could be translated in a sum of vectors in which the random vector has a constant module and changeable direction inside a range, whereas the gradient vector module 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:

The gradient vector instead of experiencing a change in its magnitude, this will be always the unit and its direction that of the greatest rate of increase of the pheromone concentration. A random direction vector with constant module will not be literally considered, but a random turning angle starting from the gradient vector direction. Attending to the previous question how do we include the saturation effect in the resulting moth shift vector?, here the answer: the dependence on the moth saturation level (interrelated with the pheromone concentration in the field) will state in the random turning angle.

NetLogo is an agent-based programming language and integrated modeling environment. NetLogo is free and open source software, under a GPL license.

NetLogo is an agent-based programming language and integrated modeling environment. NetLogo is free and open source software, under a GPL license.