- Goals
- Enhanced Nicholson-Bailey Model
- Matlab Implementation of the Simplest IDE
- An Integro-Difference Approach to Nicholson-Bailey
- Spatially Structured Environments

- Learn how MATLAB implements movies and loops.
- Couple dispersal models with discrete-time population dynamics.
- Learn some ways to implement environmental heterogeneity.
- Begin to see some of the unique behaviors associated with spatio-temporal dynamics.

For purposes of demonstration we will look at version of the
Nicholson-Bailey discrete model for the interaction between host and
parasite species. A nice, detailed account is presented in *Mathematical Models in Biology*, pp. 79-89, by Leah Edelstein-Keshet.
The purely temporal model tracks the population of the host species, ,
and the parasite species, as a function of discrete time, . The
Enhanced Nicholson-Bailey model (ENB) maps current population levels to
future population levels according to

The right-hand-side of these equations can be viewed as a product of probabilities of encounter related to the `Law' of Mass Action, that is, the probability of an encounter between host and parasite is proportional to . The probability of a host encountering parasites twice is and the probability of encountering parasites is . If is the per parasite probability of parasitization, the probability of a host being parasitizing at least once is

The various fractions occur because the

There are three fixed points for ENB:
and a mixed
population which must be calculated using a nonlinear root
finder on the equations

This mixed population is stable for an intermediate range of fecundity () and parasite effectiveness (). When is sufficiently small the mixed population is unstable and the parasites go extinct ( ). When is sufficiently large the mixed population is also unstable, but in this case there can be stable limit cycles, oscillations leading to extinction, or period-doubling bifurcations leading to chaos.To see discrete population dynamics in MATLAB you will need to implement a

% % ENB -- a small script for implementing the Enhanced Nicholson Bailey Model % n=1; a=4; r=1.75; k=1; % set parameters for ENB ngens=50; % set number of generations h=zeros(1,ngens); p=h; % initialize population vectors p(1)=.1; h(1)=.9; % set initial conditions for t=1:ngens-1 % iterate ENB for next ngens-1 generations h(t+1)=h(t)*exp( r*(1 - h(t)/k) - a*p(t) ); p(t+1)=n*h(t)*( 1 - exp(- a*p(t) )); end % of iteration over generations plot(h,'b'),hold on, plot(p,'r'),hold off legend('host','parasite')

Notice the `for` and `end` statements above, which define a loop in
the MATLAB script.
Try running this program for increasingly large values and see how the dynamics changes.

__EXERCISE 13__: The (un-enhanced) Nicholson-Bailey model for host-parasite
interactions is

Earlier (lab 1) we examined the stability of the single non-zero fixed point for this model and showed that it is always unstable. Write a MATLAB M-file which will allow you to set initial conditions exterior to the program and then iterate for some number of generations. (

Probably the simplest possible integro-difference equation is one which incorporates a single step of dispersal followed by a single step of linear, Malthusian growth at a per-capita rate of :

The order of these two steps (reproduction/dispersal or dispersal/reproduction) is not important, mathematically, although it can have significant biological repercussions (see Anderson (1991) for a discussion in the context of plants). One may think of this as a model for an insect whose adult females disperse () and then lay eggs with an overall fecundity or (including number of eggs laid per female, fraction of eggs which hatch, fraction of larvae which survive to adulthood). The following MATLAB script executes this scenario for simple drift and diffuse dispersal over several time steps.

% % MALTHUS % % matlab code to evaluate dispersal due to convolution % of a spatial dispersal probability with a population density % function. % % 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) % c - center of spatial dispersal % D - diffusion constant (per step) % r - per-capita growth constant % np=128; xl=5; dx=2*xl/np; c=.4; D=.25; r=1.15; nsteps=10; % % discretize space, define an initial population distribution % x=-xl+[0:np-1]*dx; p=(abs(x+3)<=.5); p0=p; % keep track of initial population for comparison % % define a dispersal kernel % K=exp(-(x-c).^2/(4*D))/sqrt(4*pi*D); % % normalize K to make it have integral one % 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. The fftshift accounts for % using an interval of (-xl,xl) as opposed to (0,2*xl). % fK=fft(fftshift(K))*dx; % % begin iterating the model % for j=1:nsteps, % % calculate the fft of p; perform the convolution % fp=fft(p); fg=fK.*fp; % % fg now contains the fourier transform of the convolution; % invert it, multiply by post-dispersal reproduction, % and take a look. % g=r.*abs(ifft(fg)); plot(x,p0,'r',x,g,'b'), axis([-xl xl 0 1.6]) pause(.5) % pause for .5 sec to see the plot % % update p and move to the next time step % p=g; end % of the iterationsThere are some features to note about this script. The

hold on, plot(x,p0,'r',x,g,'b'), axis([-xl xl 0 1.6]), hold offThis should allow you to see all the iterations on the same plot, with one additional per each `click' of the space bar.

__EXERCISE 14__: Imagine now that the net fecundity, , depends on space - that is,
eggs hatch and survive better in some regions or others. Modify the program
above so that is 1.3 for and elsewhere. (*Hint*:
the MATLAB command `ones(1,np)` generates a row-vector of ones of length
`np` and the logical statement `(abs(x) <= 1.5)` generates a vector
which has elements 0 where and 1 where ). Will the
population persist in this case? Now try the same scenario with with mean
dispersal distance zero ().

The idea for implementing an IDE approach to known discrete dynamics is to imagine that the behavior of the organisms can be represented in two stages:

**Reproduction and lifecycle dynamics**. This is the portion of the populations' behavior that is already captures in the ENB, equations (1,2). For this particular model, we think of parasitized hosts as immobilized, and future hosts as eggs (also immobilized). The major change is that now and as being population densities depending on space, and .**Dispersal**. Cuddly fuzzy baby hosts and parasites disperse after birth, each according to species-specific dispersal probabilities, and respectively. The basic form of the dispersals for each species will be a form of dispersal kernel derived from the heat equation (for random diffusion in a plane):

Here is the dispersal time of each species and is a species parameter describing the `diffusivity,' or rate of dispersal, of each population.

- Calculate the number of wrinkly baby hosts and parasites based on
the last population values:

- Disperse the toddling baby organisms using convolutions as described
above:

Below appears psuedo-code which implements the above model in MATLAB. For interesting results parameter values are chosen where oscillatory instabilities and chaos occur.

% % NBSPACE % % matlab code to iterate an enhanced Nicholson-Bailey model % for host and parasite, with each capable of diffusive movement. % % Intermediate variables (*) are created using a discrete % enhanced N-B step: % % h* = h(t) exp[ r (1 - h(t)/k) - a p(t) ] % p* = n h(t) ( 1 - exp[- a p(t) ] ) % % This population is then allowed to disperse using a spatial % convolution, * below, (like a long step of a heat equation): % % h(t+1) = p_h(x,y) * h* % p(t+1) = p_p(x,y) * p* % % p_h and p_p can be any probability functions representing the % dispersal of the two species; I have used Gaussians below. % % Parameters appearing (and their interpretation): % % r reproduction rate of host species % k carrying capacity of host species % a encounter rate or effectiveness of parasite on host % n parasites produced by a successful infestation % % Parameters used for the code: % % np number of points in x, y directions (a power of 2) % xl length of domain in x and y directions % dx grid spacing % dt a time-step length % mu followed by h and p -- diffusion parameter for each % species % ngens number of steps to run the code (number of % generations % % set parameters and spatial grids % np=64; mup=.02; muh=.02; dt=.5; xl=10; dx=2*xl/np; ngens=20; x=linspace(-xl,xl-dx,np); y=x; [X,Y]=meshgrid(x,y); % % set up spatial parameters % n=ones(np); a=4*ones(np); r=1.75*ones(np); k=ones(np); % % Set up stuff for movie % % M=moviein(ngens); % If you want movies % % Set up initial conditions % p0=.5*(1+cos(.5*pi*X/xl+pi/4*rand(np)).*sin(pi*Y/xl-pi/2*rand(np))); h0=.5*rand(np).*(1+cos(4*pi*sqrt(X.^2+Y.^2)/xl)); % % Define movement kernels for host and parasite % hker=exp(-(X.^2+Y.^2)/(4*muh*dt))/(4*pi*dt*muh)*dx^2; Fhker=fft2(hker); % FFT is taken because we will use it often, % two factors of dx; one for each dimension pker=exp(-(X.^2+Y.^2)/(4*mup*dt))/(4*pi*dt*mup)*dx^2; Fpker=fft2(pker); % % Now we are in a position to iterate for % a number of generations equal to ngens % for j=1:ngens % % Advance the life cycles using adjusted N-B model % hn=h0.*exp(r.*(1-h0./k)-a.*p0); pn=n.*h0.*(1-exp(-a.*p0)); % % Do the dispersal step % fhn=fft2(hn); % First take fft of each species fpn=fft2(pn); % First take fft of each species % % Now do the convolutions. % The fftshift serves to center the probability functions. % h0=real(fftshift(ifft2(Fhker.*fhn))); p0=real(fftshift(ifft2(Fpker.*fpn))); % % Now plot the solution % pcolor(X,Y,h0-p0), shading interp, colormap(jet), axis square, colorbar % M(:,j) = getframe; % If you want movies end % % Here ends the iteration % % % To play the movie, use this command: % movie(M,2,8)

Now, if you want to run this M-file as it stands you should see basically a complicated spatial field at the end. This field is coded so that the difference between relative abundances of host and parasite can be seen; places with more abundant hosts will appear as red, more abundant parasties as blue. To see a movie you must un-comment the

% If you want movieslines in the script. These movies can be memory-intensive, so don't get freaked out if your computer suddenly starts whirring, clicking and smoking! Also, movies are created in MATLAB basically by saving screen-grabs of the (normally computationally intensive) graphics. This allows the frames of the movie to be `stacked' together and viewed rapidly, but it also means that you can not change the size of movie when it is being viewed, and if you pop up any windows which overlap with the MATLAB graphics window while you are producing a movie, MATLAB will grab whatever you popped up and keep it in the movie! (I had a particularly bad moment discovering this making a movie for a conference. While producing the movie I happened to go lingerie shopping for my wife at the Victoria's Secret web site.....) The command

movie(M,2,8)will play the movie,

__EXERCISE 15__: Copy the above MATLAB M-file and modify it to implement an IDE for the
normal, vanilla Nicholson Bailey model. Run the model with initial conditions
which are:

- Spatially homogenous (constant).
- Spatially random (remember
`rand`!). - On a spectrum of structures (try taking a constant +
sinusoidally varying perturbations and then change the relative magnitudes -
as in
`k=1+.9*cos(4*pi*X/xl).*cos(4*pi*Y/xl)`).

__EXERCISE 16__: Create a MATLAB M-file which will simulate the spatial spread of a single population satisfying (when space is not included) a logistic model:

with random dispersal. For a particular (small) diffusion parameter, investigate the effect of increasing the intrinsic growth rate () of the population. What do you observe? Is this reasonable? What happens when enters the chaotic regime?

One of the nice things about IDE models implemented in MATLAB is that they can
accomodate spatial structure easily. Returning to the ENB model for a moment,
consider the ways in which the spatial environment might influence parameters
in the model. It is possible that motility of organisms can vary with the
environment - an example of this is provided above in the section on
victim-taxis. Neglecting this, the environment can have a plethora of
life-history effects, perhaps the largest of which is alteration in the
`carrying capacity,' , of the environment for the host. Returning to the
MATLAB code for the ENB given above, notice that when parameters are defined
they are defined as *matrices* of the same size as all spatial variables.
When the parameters are used in the life-history step of the IDE, the MATLAB
implementation accounts for element-by-element matrix multiplication.
Consequently, if we wanted to simpulate the effect of a crop on the host, we
could insert

k=.1*ones(np)+(abs(X)<=.25)+(abs(X)>=1.75 & abs(X)<=2.25)... +(abs(X)>=3.75 & abs(X)<=4.25)... +(abs(X)>=5.75 & abs(X)<=6.25)... +(abs(X)>=7.75 & abs(X)<=8.25)... +(abs(X)>=9.75 & abs(X)<=10.25);in place of the existing line

k=ones(np)(Notice the MATLAB form of `continuation' for a long line: the ellipsis `...') This will have the effect of putting `stripes' of high carrying-capacity crops in the environment.

__EXERCISE 17__: Implement the crop-structured spatial heterogeneity indicated above.
Think of a way to introduce parameters so that the width of the crop stripes
and their relative `value' to the host are easy to adjust. Try some
simulations - are there `crops' and spatial structures which are more or less
vulnerable to infestation by the host even in the presence of the parasite?

__EXERCISE 18__: Suppose that the cropping structure also influences the `efficiency,'
, of the parasite. How should this be included, and what are its effects
on the dynamics?

Although they may be tedious, artificial spatial structures are much easier to
introduce than `natural' ones, just because a semi-structured `random'
environment is very hard to put together off the top of your head. Of course,
in principle one could take field measurements of a spatial arena and then
scan these into density plots, but probably that is beyond the scope of our
current lab experience. One alternative was suggested to me by Professor Jon
Allen (Department of Entomology and Nematology, University of Florida,
Gainesville). The following MATLAB script, named `rockies.m' generates a
randomized but spatially correlated structure by choosing random phases and
angles in Fourier spectrum, but using a user-specified power law for
decreasing energy in the power spectrum. (If this doesn't make sense to you,
don't sweat it; the program will work anyway). The MATLAB script requires a
`fractal dimension', `H`, which increases the
spatial correlation (think regularity) as the `dimension' increases. The
code appears below, but is available for download in the class directory. The
script will produce a surface plot of the field `B`, which is normalized
so that its peak value is one, and is also then available to be used for
simulating
resource heterogeneity. The `power-law'
decrease in the spectral energy of waves, `p`, is defined as `p=(H+1)/2`
because when the exponent is 1/2 the distribution is theoretically
indistinguishable from `white' noise. The most interesting structures appear
for `H` between .5 and 1.5.

% % ROCKIES % % This cryptic matlab script generates a structured but random spatial % environment of size NxN. The user inputs the following: % N - number of grid points in both x,y directions % H - fractal dimension (the bigger this is the smoother the landscape) % N=64; H=1; p=(H+1)/2; [I,J]=meshgrid([-N/2+1:N/2], [-N/2+1:N/2]); % wave modes randn('seed',0); phases=2*pi*rand(N); % random phases amp=randn(N).*(I.^2+J.^2).^(-p); % normal amplitudes about power law amp(N/2,N/2)=0; % mean amp of zero fA=amp.*exp(i*phases); % field in tranform space A=ifft2(fA); % invert transform B=abs(A); % make field real B=B/max(max(B)); % normalize surf(B), colormap(jet), shading interp title(['Log(PS) = - ',num2str(p),' k']) view(-37.5,60)

__EXERCISE 19__: Download the program `rockies.m` and run it for a few choices of the
`fractal dimension.' Describe the relationship between increasing this
parameter and the output field.

__EXERCISE 20__: Use the `B` field output by the
`rockies` program to generate random
landscapes of carrying capacity for the ENB model. One way to do this would
be simply to execute the rockies script and then state

k=2*B;in the ENB program. Modify the graphic output to include contours of the resource with the density plots of the relative abundance of victims/parasites. What dynamics behaviors seem to correlate with what landscape features? What consequences do you think this might have for Nicholson's hypothesis for the persistence of parasitized populations in space?

*Wow, Good Job! You now know secret and official things about the
computational implementation of IDE using MATLAB !*

2001-04-06