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.
The general format of an IDE for a population is
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 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 . 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 ). 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?
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,
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 (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 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 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 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 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
% % 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!