Subsections

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

Lab 3 - Implementing an Integro-Difference Model

Goals

In this lab we will learn the philosophy and some of the practical MATLAB pitfalls associated with implementing an Integro-Difference Equation (IDE) model. An IDE is one natural way to add spatial dynamics to a temporally discrete model. In this lab you will:

• 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.

Enhanced Nicholson-Bailey Model

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, $H_t$, and the parasite species, $P_t$ as a function of discrete time, $t$. The Enhanced Nicholson-Bailey model (ENB) maps current population levels to future population levels according to

$$H_{t+1} = f(H_t, P_t) = H_t \exp \left[ r \left(1 - \frac{H_t}{K} \right) - a P_t \right],$$ (1)

$$P_{t+1} = g(H_t, P_t) = n H_t \left( 1 - \exp[-aP_t] \right) .$$ (2)

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 $H_t \cdot P_t$. The probability of a host encountering parasites twice is $H_t \cdot P_t^2,$ and the probability of encountering $n$ parasites is $H_t \cdot P_t^n$. If $a$ is the per parasite probability of parasitization, the probability of a host being parasitizing at least once is

$$\underbrace{a H_t P_t}_{\mbox{\tiny 1 encounter}} + \underbr... ...ny 3 encounters}} + \cdots = H_t\left[ 1 - e^{-a P_t}\right] .$$

The various fractions $\frac1{n!}$ occur because the order of the encounters does not matter. For example, there are $6= 3!$ ways for three different parasites to encounter a particular host, but if the order of these encounters does not matter we need to factor out that $3!$ to get the probability of three order-independent encounters. Consequently, the right-hand-side of $\ref{psiteeq}$ can be thought of as the expected number of hosts which have been parasitized at least once times the number of fuzzy baby parasites ($n$) produced by each successfully infected host. The product $H_t \exp[-aP_t]$ on the RHS of equation (1) is the expected number of non-parasitized hosts, which then reproduce with fecundity based on the Ricker map, $\exp\left[ r \left(1 - \frac{H_t}{K} \right) \right]$.

There are three fixed points for ENB: $(H,P) = (0,0), (K,0)$ and a mixed population $(H^*,P^*)$ which must be calculated using a nonlinear root finder on the equations

$$H^*=f(H^*,P^*) \quad \mbox{and} \quad P^*=g(H^*,P^*).$$

This mixed population is stable for an intermediate range of fecundity ($r$) and parasite effectiveness ($a$). When $a$ is sufficiently small the mixed population is unstable and the parasites go extinct ( $H_t \rightarrow K$). When $a$ 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 for loop. For example, to see behavior in the ENB you can implement the following M-file:

     % 
     %   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 $r$ values and see how the dynamics changes.

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

$$\begin{eqnarray*} H_{t+1} &= &\lambda H_t e^{-aP_t}, \\ P_{t+1} &=& c H_t \left(1-e^{-aP_t}\right). \end{eqnarray*}$$

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. (Hint: you may wish to modify the RHS of these equations to remove the possibility of negative values!) Additionally, plot the populations against one another (phase space) as well as against time (state space). The parameters $a=0.068, c=1, \lambda=2$ correspond to the interaction between a greenhouse whitefly and its chalcid parasitoid. Now you have visualized the well-known instability of what seems like a pretty reasonable model. Nicholson (1957) proposed that the instability and subsequent death of both species would be moderated in nature by continual re-invasion of patches.

Matlab Implementation of the Simplest IDE

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 $r$:

$$\begin{eqnarray*} P_{n+\frac12} &=& K*P_n \\ P_{n+1} &=& r P_{n+\frac12} . \end{eqnarray*}$$

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 ($P_{n+\frac12}$) and then lay eggs with an overall fecundity or $r$ (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 iterations

There are some features to note about this script. The pause(.5) command causes the program to stop for half a second and plot. Try commenting out this line and running the program again - you will see only the end result as opposed to the intervening steps. Now try removing the parentheses, that is, replace pause(.5) with pause. This will cause the program to wait for user input before proceeding - for example, just hit the space bar in the command window and the program will advance one step. Finally, to see all the steps together, replace the plot line with

   hold on, plot(x,p0,'r',x,g,'b'), axis([-xl xl 0 1.6]), hold off

This 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, $r$, depends on space - that is, eggs hatch and survive better in some regions or others. Modify the program above so that $r$ is 1.3 for $\vert x\vert<1.5$ and $r=.8$ 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 $\vert x\vert>1.5$ and 1 where $\vert x\vert\le 1.5$). Will the population persist in this case? Now try the same scenario with with mean dispersal distance zero ($c=0$).

An Integro-Difference Approach to Nicholson-Bailey

Description of IDE for ENB

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 $H_t$ and $P_t$ as being population densities depending on space, $H_t(x,y)$ and $P_t(x,y)$.

• Dispersal. Cuddly fuzzy baby hosts and parasites disperse after birth, each according to species-specific dispersal probabilities, $K_H(x,y)$ and $K_P(x,y)$ 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):
  
  $$K(x,y,T) = \frac{1}{4 \pi \mu T} \exp\left[ -\frac{x^2 + y^2}{4 \mu T} \right].$$
  
  
  Here $T$ is the dispersal time of each species and $\mu$ is a species parameter describing the `diffusivity,' or rate of dispersal, of each population.

The mathematical blueprint for the IDE approach to ENB is then

1. Calculate the number of wrinkly baby hosts and parasites based on the last population values:
   $$\begin{eqnarray*} H_{t+1/2}(x,y) = & f(H_t, P_t) &= H_t(x,y) \exp \left[ r \le... ...H_t, P_t) &= n H_t(x,y) \left( 1 - \exp[-aP_t(x,y)] \right) . \end{eqnarray*}$$
   
   
   
   
   

2. Disperse the toddling baby organisms using convolutions as described above:
   $$\begin{eqnarray*} H_{t+1}(x,y) &= & K_H * H_{t+1/2}(x,y), \\ P_{t+1}(x,y) &= & K_P * P_{t+1/2}(x,y) . \end{eqnarray*}$$
   
   
   
   
   

Notice that the dispersal kernels may be different for host and parasite! By iterating these rules for many $t$ we can simulate many generations of reproduction, interaction, and dispersal.

Matlab Implementation

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 movies

lines 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, M, 2 times at a speed of 8 frames per second. The advantage of having a movie is that you can play it again and again without re-doing the calculation. You can also save a movie to a transportable AVI file.

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:

1. Spatially homogenous (constant).
2. Spatially random (remember rand!).
3. 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)).

Can you verify Nicholson's hypothesis that space allows persistence of both species?

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:

$$P_{n+1} = P_n + r P_n (1 - P_n)$$

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

Spatially Structured Environments

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,' $K$, 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,' $a$, 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 !