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

- Learn how to use even and odd reflections to implement lethal ()
and solid (
) boundary conditions for IDE.
- Implement Powell et. al.'s 1998 approach for modelling chemotactic
movement of a population.
- Spend time working out some of the computational details of the student
case studies.

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.

The general format of an IDE for a population is

Here is the functional response connecting the dispersal stage organism (with density ) to the previous generation and 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 at some
. 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, , 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 ), the effect of can be simulated by
creating an *anti-*population of dispersers on the opposite side of ,
which when it disperses back toward the positive (real) population will
exactly cancel out at . Since the functions we are looking at must be
periodic on the entire (reflected) interval, we will also have a lethal
boundary at 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

__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?

A reflecting boundary (that is, a boundary condition
) 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.

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 and a
lethal (absorbing) boundary at . 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 dispersalNotice 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.

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 and ). Let's implement reflecting
boundaries for a `ring-random' or `ripple' dispersal kernel,

which models propagules leaving the place of their origin at speed and travelling and dispersing (with diffusivity ) for a time 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 must be evaluated numerically for normalization. The MATLAB script below implements `ripple' dispersal with reflecting boundaries at and and lethal boundaries at and . First, we set up both the grid of interest,

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 convolutionThis 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 direction for the reflecting boundaries:

py=[p0; p0(np:-1:1,:) ];The semi-colon above places one matrix above the other, and the

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 squareThe 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 and indices from

__EXERCISE 23__: To make sure that you understand the reflection process, edit the script
above so that it implements a lethal boundary at and a
reflecting boundary at and . How many reflections would you need
to implement a reflecting boundary at *only* , with lethal boundaries
on the other three?

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 , fermenting and
releasing ethanol () at a rate . The steady-state ethanol
distribution then satisfies

where 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

Here is a saturation parameter - basically the level at which the sensory apparatus of the

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

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

The constant must be chosen to normalize 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 iterationsOne can think of this implementation as a pre-multiplication of the dispersing population, , with a dispersal `bias': individuals further away from the source of the ethanol (therefore sensing low ) are more strongly biased to random dispersal. The biased population () is dispersed at random, and then de-biased (by multiplication of ). 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 - five steps with is not necessarily equivalent to two steps with .

__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 () 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!*

2001-04-06