next up previous contents
Next: Bibliography Up: Spatio-Temporal Models in Ecology; Previous: Lab 3 - Implementing

Subsections


Lab 4 - Boundaries, Chemotaxis and Student Case Studies

Goals

In this lab we will finish off some of the complicated details which may be necessary for practical application of IDE. The first of these is the implementation of boundaries, in particular `lethal' (or `absorbing') boundaries and `solid' (or `reflecting') boundaries. The second is accounting for semi-random movement in response to external stimuli, in particular chemo-taxis. Thus, we will

An overall goal for this class is to get students to `work without a net,' that is, to create IDE simulations of ecological circumstances which may shed light on research questions. To encourage this, we will have lab time for teams to begin implementing an IDE approach, with our fearless instructor around to help with difficulties. Participants should work in teams of 2 or more to discuss the models which will be implemented, parameter regimes to be studied, and experimental protocols to be used. Then each student should build their own simulator and implement their portion of the experimental protocol, and get back together with the rest of the group to discuss results. I have suggested some group projects, but feel free to work on something closer to your own research or interests. Each group should have something to present to the rest of the class on the final day. You may have to work mainly with one spatial dimension in order to have enough time to experiment freely. As we have seen in the rest of the class, if you can do it in one dimension with MATLAB then you can do it in two with only slight modification.

Reflecting and Absorbing Boundary Conditions

The Wall of Doom

The general format of an IDE for a population $P$ is

\begin{eqnarray*}
P_{t+1/2} &=& f(P_t), \\
P_{t+1} &=& K  *  P_{t+1/2}.
\end{eqnarray*}



Here $f$ is the functional response connecting the dispersal stage organism (with density $P_{t+1/2}$) to the previous generation and $K$ is the spatial dispersal probability for an initially localized individual. Spatial boundaries only have an impact in the second step during which dispersal occurs.

A lethal boundary condition is written mathematically as $P=0$ at some $x=L$. This would seem to be quite a problem for our FFT-based approach, which assumes at a basic level that functions going into the FFT are periodic. However, we can get around this by reflecting the data. If the dispersal kernel, $K$, is symmetric around the origin (so that the probability of moving to the left is exactly the same as that of moving to the right, or $K(-x) = K(x)$), the effect of $P(x=L) =0$ can be simulated by creating an anti-population of dispersers on the opposite side of $x=L$, which when it disperses back toward the positive (real) population will exactly cancel out at $x=L$. Since the functions we are looking at must be periodic on the entire (reflected) interval, we will also have a lethal boundary at $x=0$ because of this reflection. Below is a script which executes this implementation of boundary conditions.

     %
     %       BCs
     %
     %       matlab code to illustrate the effect of  boundary conditions 
     %       and their implementation
     %       
     %       boundary conditions will be implemented  at x=0, x=xl
     %       by reflecting the domain
     %
     %       p(x)    -       population density in x
     %       K(x)    -       probability of individual dispersal to location
     %                          x starting at location 0.
     %       xl	     -       maximum x coordinate (also - minimum)
     %       dx	     -       spatial step  delta x = 2*xl/np
     %       np	     -       number of points  in discretization for x
     %                         (should be a power of 2 for  efficiency of FFT)
     %       D	     -       diffusion constant (per step)
     %       
            
     np=128; xl=5; dx=2*xl/np; c=0; D=.25; 
     
     %
     %       discretize space, define an initial population distribution
     %
     
     x=-xl+[0:np-1]*dx;
     p=(abs(x-1.5)<=1.); p0=p;  %  keep track of initial population 
     plot(x,p0,'r')
     
     %
     %       define and normalize a dispersal kernel
     %
     
     K=exp(-(x).^2/(4*D))/sqrt(4*pi*D);
     K=K/(dx*sum(K));
     
     %
     %       calculate the fft of K, multiplying by dx to account
     %       for the additional factor of np and converting from a
     %       interval length of 1 to 2*xl.  
     %
     
     fK=fft(K)*dx;
     
     %
     %       skew-reflect the data to implement wall of doom
     %          pa contains the actual population before dispersal
     %          pr contains the reflected populuation before dispersal
     %
        
     pa=p;
     pr=-p(np:-1:1);  %  start at last index, go  backwards for reflection
     
     %
     %       calculate the fft of the pa, pr; perform the convolutions
     %
     
     fpa=fft(pa);     fpr=fft(pr);
     fga=fK.*fpa;     fgr=fK.*fpr;
       
     %
     %       fgr and fga now contain FFT of the convolution of reflected 
     %       and actual fields, respectively; now invert
     %

     gr=real(fftshift(ifft(fgr)));
     ga=real(fftshift(ifft(fga)));
     g=ga+gr;

     %
     %       update p, undoing the reflection (that is, take only the data
     %       from x=0 to xl, with zeros elsewhere
     %
     p=g.*(x>=0);
     
     hold on, plot(x,p,'m',x,gr,'y',x,ga,'b',x,g,'g--'),  hold off
     axis([-xl xl -1 1])
Make particular notice of the statement
     pr=-p(np:-1:1);
which implements the reflection of the data. In MATLAB the statement np:-1:1 will generate numbers starting at np and stepping backwards (by steps of -1) to 1. Placing these indices into p, as above, reflects the field (negatively) around $x=0$. The program then disperses both the original and reflected fields, re-sums to get the field which implements the boundary conditions (g above), and then recovers the actual field (by setting the field back to zero for negative $x$). The results are plotted with the reflected field in yellow, the unreflected field in blue, the summed field in green dashes, and the final result (the population after dispersal, taking into account the lethal boundaries) in magenta.

EXERCISE 21: Use for ... end statements and pause to continue the dispersal through several time steps. When you are confident that you have done this properly, extend the program to account for reproduction after dispersal (that is, replace p=g.*(x>=0) with p=r*g.*(x>=0), where r is a parameter you will set to be larger than 1). What do you observe after iterating the combination of dispersal, lethal boundaries, and reproduction for several generations?

Reflecting Boundaries

A reflecting boundary (that is, a boundary condition $\frac{\partial P}{\partial x} = 0$) is implemented the same way as a lethal boundary, except that the reflection is done positively instead of negatively.

EXERCISE 22: To implement solid boundaries, in the script you have written above replace the line

     pr=-p(np:-1:1);
with the line
     pr=p(np:-1:1);
Run the script and see what happens. If you iterate dispersal several times the final result should approach a constant, non-zero average value. Why is this? Remember that you have a solid wall on each end of the domain. The effect of a long period of random motion, since no individuals are lost, is to eventually create a uniform distribution.

Mixed Boundaries

It might happen that one boundary is lethal while the other is solid, in which case one needs to do a double reflection. Let's suppose we want to investigate a situation with a solid (reflecting) boundary at $x=0$ and a lethal (absorbing) boundary at $x=L$. A MATLAB script implementing mixed boundary conditions (with drift and drop dispersal) appears below.

     np=64;
     xl=5;
     dx=xl/np;
     x=dx*[0:np-1];            % the physical space of interest
     x2=dx*[-2*np:2*np-1];     % the virtual space on which reflections occur

     D=.4;                     % per-step dispersal distance

     K=.5/D*exp(-abs(x2)/D);   % double-exponential dispersal kernel
     fK=dx*fft(K);

     p0=(abs(x-2.5)<=1);       % initial population

     for istp=1:10,            % begin several dispersal steps

         pr=[p0, -p0(np:-1:1)];    % negative reflection to right for
                                   %    lethal boundary at x=xl

         pl=[pr(2*np:-1:1), pr];   % positive  reflection to left for
                                   %    reflecting boundary at x=0

         fpl=fft(pl);
         p2=real(fftshift(ifft(fK.*fpl)));

         pn=p2(2*np+1:3*np);       % use only that data corresponding to the
                                   %     region of interest

         hold  on, plot(x,p0,'r',x,pn,'g'), axis([0 xl 0 1]), hold off
         pause(.1)
         p0=pn;                    % update the population field

     end                       % of dispersal
Notice that implementing this kind of mixed boundaries requires vectors twice as long as implementing boundaries of the same type. This is not such a big problem in one dimensional simulations, but it can become prohibitive in multiple dimensions.

Boundary Conditions in Two Dimensions

In principle, lethal and reflecting boundaries are implemented the same way in two dimensions as in one, and to implement boundaries the dispesal kernel must still be symmetric. Now, however, there are two sets of reflections to implement (boundaries in both $x$ and $y$). Let's implement reflecting boundaries for a `ring-random' or `ripple' dispersal kernel,

\begin{displaymath}
K\left(r = \sqrt{x^2+y^2} \right) = C e^{-\frac{(r-v t)^2}{4 D t}},
\end{displaymath}

which models propagules leaving the place of their origin at speed $v$ and travelling and dispersing (with diffusivity $D$) for a time $t$ before settling. This kernel (also called the `ripple') was used by Brewster et. al. (1997) to investigate the dispersal of whiteflies in the Imperial Valley of California. The constant $C$ must be evaluated numerically for normalization. The MATLAB script below implements `ripple' dispersal with reflecting boundaries at $y=0$ and $y=L$ and lethal boundaries at $x=0$ and $x=L$. First, we set up both the grid of interest, [X,Y], and the grid required for the reflections, [X2,Y2].
     t=2; D=.1; v=2;                               % set parameters

     xl=10; yl=10; np=64; dx=2*xl/np; dy=2*yl/np;  % parameters for the
                                                   % x,y coordinates

     x=linspace(0,xl-dx,np); x2=linspace(-xl,xl-dx,2*np);
     y=linspace(0,yl-dy,np); y2=linspace(-yl,yl-dy,2*np);

     [X,Y]=meshgrid(x,y);                          % grid of interest
     [X2,Y2]=meshgrid(x2,y2);                      % reflection grid

     K2D=exp(-(sqrt(X.^2+Y.^2)-v*t)^2./(4*D*t));   % ripple dispersal
     K2D=K2D/(dx*dy*trapz(trapz(K2D)));            % normalize

     fK=fft2(K2D)*dx*dy;                           % FFT for convolution
This script sets up the grids and the ripple dispersal kernel. Now we can define an initial population (in a circle centered in the domain)
     p0=(((X-xl/2).^2+(Y-yl/2).^2)<=4);
and perform reflections preparatory to doing the dipsersal. First we must do an even reflection in the $y$ direction for the reflecting boundaries:
     py=[p0; p0(np:-1:1,:) ];
The semi-colon above places one matrix above the other, and the (np:-1:1,:) takes rows (first index) in reverse order, for all columns (the colon in the second argument). This matrix must now be negatively reflected in the $x$ direction to implement the lethal boundary
     px=[-py(:, np:-1:1) py];
Now we can just do the regular dispersal on the reflected grid:
     fp2=fft2(px);
     p2=real(fftshift(ifft2(fp2.*fK)));
     pt=p2(np+1:2*np,np+1:2*np);
     pcolor(X,Y,pt), shading interp, colormap hot, axis square
The only particularly difficult thing here is knowing which part of the reflected and dispersed data to take to get the population of interest. Given the way that we did the reflection it is the `upper-right-corner' of the data, or $x$ and $y$ indices from np+1 to 2*np.

EXERCISE 23: To make sure that you understand the reflection process, edit the script above so that it implements a lethal boundary at $y=0, L$ and a reflecting boundary at $x=0$ and $x=L$. How many reflections would you need to implement a reflecting boundary at only $x=0$, with lethal boundaries on the other three?

A Model for Chemo-Taxis

The last thing that might be of interest is implementing chemo-taxis. Many (if not most) insect species find mates, hosts, or prey using chemical clues. Bark beetles find hosts and initialize a `mass-attack' using a pheromone feedback system (Powell et. al. 1998), natural populations of Drosophila find rotting fruit to lay eggs on by following the odor of ethanol, and a variety of insects, including ladybird beetles (Coccinella) respond to the odor of sugar-water. Powell and co-workers (1998) developed an IDE approach to emulate the chemotactic dispersal process, which is discussed here in one dimension in the context of Drosophila response to ethanol.

Consider a resource (apples) distributed with density $R$, fermenting and releasing ethanol ($E$) at a rate $\delta$. The steady-state ethanol distribution then satisfies

\begin{displaymath}
E = K_E \star \left( \delta R \right), \quad K_E = \frac1{2L} \exp\left[
-\frac{\vert x\vert}{L}
\right],
\end{displaymath}

where $L$ is the mean dispersal distance of the ethanol. A model for the `sensory index', or degree of saturation of the sensory apparatus of the flies, is

\begin{displaymath}
F(E) = \frac{E}{E_0+E}
\end{displaymath}

Here $E_0$ is a saturation parameter - basically the level at which the sensory apparatus of the Drosohpila is halfway to being completely saturated with ethanol. A model for the population repsonse is

\begin{displaymath}
\frac{\partial}{\partial t} P = - \frac{\partial}{\partial x...
...\partial P}{\partial
x}}_{\mbox{\tiny random motion}}
\right].
\end{displaymath}

The parameters $\nu$ and $\mu$ paramerize the relative strengths of the tactic response and random dispersal, respectively. Let

\begin{displaymath}
K_R = \frac1{\sqrt{4 \pi \mu t}} \exp\left[
-\frac{x^2}{4 \mu t}
\right]
\end{displaymath}

be the random dispersal kernel for flies (when $\nu=0$ or $E=0$), implemented over some time interval, $t$. An approximate solution (details of how approximate are discussed in Powell et. al., 1998) is given by

\begin{eqnarray*}
\Pi(x,t) &=& K_R \star \left\{ \exp\left[ -\frac{\nu}{\mu} F(E...
...
P(x,t) &=& C \exp\left[ \frac{\nu}{\mu} F(E) \right] \Pi(x,t).
\end{eqnarray*}



The constant $C$ must be chosen to normalize $P(x,t)$ because this method is not guaranteed to preserve the number of individuals in the dispersing population. Following is a MATLAB code which implements the chemotactic procedure over several iterations.

     %
     %       CHEMOTAX
     %
     %       matlab code to emulate a chemotactic process
     %
     %       p(x)    -       population density in x (drosophila)
     %       e(x)    -       ethanol concentration
     %       r(x)    -       resource distribution
     %       KR(x)   -       probability of individual dispersal
     %       KE(x)   -       dispersal kernel of ethanol
     %       xl	-       maximum x coordinate (also - minimum)
     %       dx	-       spatial step  delta x = 2*xl/np
     %       np	-       number of points  in discretization for x
     %       mu	-       diffusion constant (per step)
     %       nu      -       chemotactic constant
     %       a       -       dispersal distance of ethanol
     %       d       -       release rate of ethanol from resource
     %       E0      -       saturation constant for sensory index
     %       
       
     np=128; xl=8; dx=2*xl/np; 
     t=1; mu=1; nu=10;
     a=.5; d=2; E0=.5;
     nsteps=15;

     %
     %       discretize space, define an initial population distribution
     %

     x=-xl+[0:np-1]*dx;
     p=(abs(x+3)<=.5);   p0=p;   %  keep track of initial p for comparison
     r=(abs(x-2)<=.5);           % resource density

     %
     %       define dispersal kernels and normalize
     %

     KR=exp(-(x).^2/(4*D*t))/sqrt(4*pi*D*t);
     KE=exp(-abs(x)/a);
     KR=KR/(dx*sum(KR));
     KE=KE/(dx*sum(KE));

     %
     %       calculate the fft of KR, KE, dx accounts for  convolution
     %

     fKR=fft(KR)*dx;
     fKE=fft(KE)*dx;

     %       calculate ethanol field

     E=real(fftshift(ifft(fft(d*r).*fKE)));
     plot(x,p,'b',x,p,'m--',x,E/max(E),'y',x,r,'r'), axis([-xl xl 0 1.1*max(p)])

     %       calculate sensory index

     FE=E./(E0+E);

     %       begin iterating the model 

     for j=1:nsteps,

     %       
     %       calculate the fft of pi = p.*exp(-nu/mu f(E)); 
     %       perform the convolution on the pi
     %
        peye=p.*exp(-nu/mu*FE);
        fpi=fft(peye);
        fg=fKR.*fpi;

     %
     %       fg now contains the fourier transform of the convolution;
     %       invert it, multiply by the inverse exp(nu/mu*FE)
     %

        g=exp(nu/mu*FE).*real(fftshift(ifft(fg)));
   
     %
     %       Now we must be careful to normalize!
     %
        C=trapz(p)/trapz(g); g=C*g;
   
        hold on, plot(x,g,'b'), hold off
        pause(.1)

     %       update p and move to the next time step

        p=g;
     end               % of the iterations
One can think of this implementation as a pre-multiplication of the dispersing population, $P$, with a dispersal `bias': individuals further away from the source of the ethanol (therefore sensing low $f(E)$) are more strongly biased to random dispersal. The biased population ($\Pi$) is dispersed at random, and then de-biased (by multiplication of $\exp[\nu/\mu
f]$). Provided the time steps are small enough, this approximates the chemotactic dispersal process. As it turns out, the model works just fine for large time steps. As with the prey-taxis discussed in the previous lab (Turchin-taxis) the main difficulty is that the process is not linear in $t$- five steps with $t=2$ is not necessarily equivalent to two steps with $t=5$.

EXERCISE 24: Form hypotheses about the impact of the parameters on chemotactic dispersal, and test your hypotheses against model iterations. For example, for a given time step, if the ethanol disperses further (larger a) the flies should find the resources more rapidly. Or, if chemotaxis is stronger than random motion ($\nu > \mu$) then the resulting pattern of aggregation should be much more tightly centered on the resource. Change parameters to reflect your hypotheses and test if they are correct.

Well, look at you! You have gotten to the end of the formal lab material on integro-difference equations! Good Work!


next up previous contents
Next: Bibliography Up: Spatio-Temporal Models in Ecology; Previous: Lab 3 - Implementing   Contents
James Powell
2001-04-06