Subsections

• Goals
• Probability and Dispersal in One Space Dimension
  • Dispersal via Random Walks and the Diffusion Equation
  • Convolutions and Fourier Transforms
  • Numerical FFT
  • Matlab Implementation
  • Dispersal and Settling with a Localized Source
  • Ballistic Dispersal
  • Movement at Constant Speed from a Source and Settling
• Probability and Dispersal in Two Space Dimensions
  • Dispersal via Random Walks in Two Dimensions
  • Modal Dispersal in Two Dimensions
  • Turchin's Model for Prey-Taxis

Lab 2 - Dispersal

Goals

During this lab students will use MATLAB to become familiar with methodology for simulating the probabilistic dispersal of organisms using FFTs (Fast Fourier Transforms). In particular, students will:

• Visualize and understand a variety of dispersal kernels in one and two dimensions.

• Be aware of the relationship between probability kernels and special solutions to partial differential equations (PDE).

• Use MATLAB `M' files to simplify simulations.

Probability and Dispersal in One Space Dimension

Dispersal via Random Walks and the Diffusion Equation

Many of the basic probability kernels associated with population dispersal find their basis in partial differential equations. For example, the normal distribution is associated withthe diffusion equation

$$u_t = D u_{xx}.$$

The diffusion equation, among other things, models the probability ($u$) of finding an individual near any spatial locus ($x$) at some time $t$, given that the individual is moving by taking randomly chosen steps of size $\lambda$ to the left or right at a rate of one step per time interval $\tau$. The step size and time interval occur in the diffusion equation as

$$D=\frac{\lambda^2}{2 \tau}.$$

For the dispersal of a population, the important solution to the diffusion equation is the one which begins with a perfectly localized individual at the point $x=0$ (in math we know this as the Dirac delta function, $\delta(x)$), so that the initial condition is $u(x,0) = \delta(x)$. The corresponding solution, $K_D$, is called the fundamental solution and has the form

$$K_D(x,t)= \frac1{\sqrt{4 \pi D t}} \exp\left[ \frac{-x^2}{4 D t} \right].$$

This function can be interpreted as the pdf associated with the location of an individual at time $t$, moving under random walks, which was initially located at $x=0$.

Now, suppose a population of organisms is dispersed over space with an initial density of $P(x,0)$. One can think of the number of individuals located in a small interval of size $dy$ at a location $y$ as $P(y,0) dy$. These indidviduals disperse randomly according to the distribution $K_D$; thus, their probabilistic location at a later time will be

$$P(y,0) K_D(x-y,t) dy,$$

where the dispersal kernel has had its argument shifted so that it is now centered at the original locus of individuals ($x=y$). The total population after a time $t$ would then be the sum over all such infinitesimal intervals containing populations:

$$P(x,t) = \int_{-\infty}^{\infty} P(y,0) K_D(x-y,t) dy \stackrel{\mbox{\tiny def}}{=} P(x,0) * K(x,t).$$

This latter operation is called the convolution and is defined by the preceding integral.

Convolutions and Fourier Transforms

The Fourier Transform of a function defined on an interval from $-\infty$ to $\infty$ is

$$\hat{f} (k) = \frac1{\sqrt{2 \pi}} \int_{-\infty}^\infty f(x) e^{-i k x} \ dx .$$

The transformed function $\hat{f} (k)$ can be thought of as the amount of `energy' the function has stored in the `wave' $e^{ikx} = \cos(kx) + i\sin(kx)$. A convenient property of the Fourier Tranform (FT) is that the FT of the convolution is the product of the FTs,

$$\widehat{f*g} (k) = \frac{1}{2\pi} \hat{f}(k)\hat{g}(k).$$

Consequently, a slick analytic way to evaluate a convolution without actually doing the nasty integral is:

1. Calculate the FT of $f$ and $g$.

2. In wave space, evaluate the product of $\hat{f} (k)$ and $\hat{g} (k)$.

3. $\widehat{f*g} (k) = \frac{1}{2\pi} \hat{f}(k)\hat{g}(k).$

4. Invert the transform of $\widehat{f*g}$ to obtain $f*g$.

The only turd in this particular punchbowl is that FTs are not all that easy to calculate or invert. Fortunately, there is a very simple computational routine, called the Fast Fourier Transform (FFT) which will perform these operations numerically.

Numerical FFT

The numerical (or Discrete) Fourier Transform (DFT) is defined with summations instead of integrals,

$$\hat{f}_k = \sum_{j=-N/2}^{N/2} f(j \Delta x) e^{2 \pi i k j \Delta x},$$

where the function $f$ is thought of as being defined on the interval from $x=0$ to $x=1$, chopped up into $N$ chunks of length $\Delta x$. If $N$ is a power of two, (that is $N=2^m = 2,4,8,16,32,64,128,\cdots$) there is an efficient numerical algorithm for calculating the DFT, called the FFT. One curiosity of the FFT is that (depending on the implementation) it generates a factor of $N$ in the output, that is, for input $f(j \Delta x)$ it gives output $N \hat{f}_k$ instead of simply $\hat{f}_k$. Consequently, the FFT of a convolution of two functions defined on the interval $x \in [0,1]$ is

$$\widehat{f*g}_k = \frac1N \hat{f}_k \hat{g}_k .$$

Suppose now we have a population function $p$ defined on a periodic interval $x \in [-L , L]$, and a dispersal function, $f$ defined on the same interval. A blueprint for calculating the dispersed population $f*p$ using the FFT follows:

• Use FFT to calculate $\hat{f}_k$ and $\hat{p}_k$.

• Calculate $\widehat{f*p}_k = \frac{2L}{N} \hat{f}_k \hat{p}_k$. The extra factor of $2L$ allows for a change in the length of the interval from 1 to $2L$. As a note of interest, the fraction $\frac{2L}{N}$ is equal to the discrete step size in space, $\Delta x$.

• Invert the FFT to arrive at $f*p (x)$.

An added wrinkle is that, although it is theoretically natural to center the interval at $x=0$, in discrete terms this is unnatural. In general a shift has to be performed on the discrete version of $p$ so that the value of $p$ at $x=0$ is associated with the first element of a vector, as in the following diagram:

Matlab Implementation

We can use FFTs in MATLAB to calculate convolutions rapidly. Let's begin by illustrating the procedure for a population initially localized uniformly between -1 and 1, dispersing over time into a much larger domain. The FFT assumes periodicity, so if we define a space interval between -10 and 10 with 256 nodes (remember the FFT is much more efficient for computations involving $2^m$ elements), we need to take into account that the node corresponding to 10 must be left off (the FFT assumes that whatever data might be specified there is identical with the data at $x=-10$). So, to define this space:

     xlr=10; np=256; dx=2*xlr/np;
     x=dx*[0:np-1]-xlr;        % OR, equivalently
     x=linspace(-xlr,xlr-dx,np);

The independent variable is now defined. Let's also define the normal dispersal kernel and the initial population:

     p0 = (abs(x) <= 1);
     t=1; D=1; 
     K=1/sqrt(4*pi*D*t)*exp(-x.^2/(4*D*t));

Having the relevant functions defined we can now take the Fourier transform

     fp0=fft(p0);
     fK=fft(K);

Using the fact that the transform of the convolution is the product of the transforms, we can now evaluate the convolution:

     fp1=fK.*fp0;
     p1=dx*fftshift( ifft( fp1 ) );

The vector p1 now contains the dispersed population after 1 time unit, but there are some weird things in the command that need explaining. First off, the command fftshift basically chops a vector in half and interchanges the first and second halves as blocks. In this context it is necessary because the MATLAB FFT is built on the assumption that functions are defined on the interval from 0 to 1, as opposed to centered around 0. Secondly, the factor of dx. Firstly, for numerical efficiency the FFT multiplies a vector by np, and the inverse FFT divides by that factor. Since we have implemented a product of FFTs, we need to divide out a factor of np. Secondly, numerical FFTs assume that functions are defined on an interval of length 1, but our function is defined on an interval of length 2*xlr, and so the FFT must be scaled by the interval length. Together, these two conditions are equivalent to multiplying by the step size, dx=2*xlr/np, which is convenient from our perspective.

Now let's see what we have done:

     plot(x,p0,'r',x,p1,'y')

You should also get a warning about imaginary numbers being ignored. This occurs because the FFT requires complex numbers, and numerical rounding creates very small imaginary components, even for perfectly real data fields. The best way to get around this is to take the real part after convolution, which we will do in the future. Now it is possible to keep iterating the random dispersal of the population. We can get the population after the next time step using a single long, cryptic MATLAB command:

     p2=dx*real( fftshift( ifft( fft( p1 ).*fK ) ) );

Use hold on to plot this against the previous results. Does it look reasonable? Use the up-arrow ($\uparrow$) and edit the previous command, plot, and see the population after dispersal over three time units.

EXERCISE 6: Use MATLAB to define the population density function

$$P(x,0) = \left\{ \begin{array}{lr} 2 - 2\vert x\vert, & \quad -1 \le x \le 1, \\ 0 , & \vert x\vert>1, \end{array} \right.$$

on the interval from -15 to 15. In a single window plot random dispersal of this population after times of .5, 1, 2, and 4 with diffusion constant $D=2$. (Using $\uparrow$ here will simplify your life!) Use trapz to approximate the integral of the population (total number of individuals) and convince yourself that the number of individuals is being conserved.

Dispersal and Settling with a Localized Source

Another interesting dispersal model involves assuming random motion for propagules (seeds, spores, pollen, motile larvae) while they are in the medium and then settling onto a substrate at some rate $\lambda$. If propagules are initially released at the point $x=0$ the PDE modelling the density of propagules in the medium is

$$u_t = D u_{xx} - a u, \quad u(x,0)=\delta (x).$$

The deposition of these particles on the ground is tracked by the variable $v$, which satisfies

$$v_t = a v, \quad v(x,0)= 0.$$

The long-time ( $t\rightarrow \infty$) solution for $v$ is a dispersal kernel, and can be written (following Neubert et. al. 1995)

$$K_S(x) = \sqrt{\frac{a}{4 D}} \exp\left[ -\sqrt{\frac{a}{D}} \vert x\vert \right].$$

The following commands plot this function on the previously defined x interval:

     D=1; a=2;
     KS=sqrt(.25*a/D)*exp(-sqrt(a/D)*abs(x));
     plot(x,KS,'r')

Now try plotting KS on the same axes, varying the settling rate, $a$, to get an idea of how it behaves.

EXERCISE 7: For a uniform distribution of sources between -1 and 1 investigate the dispersal of propagules according to the probability distribution $K_S$ given above. On the same set of axes plot the effect of changing $a$ on the dispersal of the population.

Ballistic Dispersal

A model for ballistic dispersal of spores from ground level in random directions is given by Neubert, Kot and Lewis, 1995, by

$$K_B= \left\{ \begin{array}{lr} \frac1{\pi \sqrt{\frac{c^4}{g... ...{g}, \\ 0, & \vert x\vert> \frac{c^2}{g}. \end{array} \right.$$

Here $c$ is the speed at which spores are ejected and $g$ is the gravitational constant. Under the assumption that every spore germinates and the sporulating adults die each generation, the spread of a population of fungi is modelled by

$$P_{n+1}(x) = K_b * P_n(x),$$

where $n$ denotes the generation number and $P_0(x)$ denotes the initial distribution of the colony.

We could illustrate the spread of the colony using the commands we have already mastered in MATLAB , but this is a good place to use an `m' file, which is a MATLAB script which runs a sequence of commands. In the File menu of the MATLAB command window click on New and then M-file. An editor/debugger window will pop up, in which you can type in a series of commands, save it as an M-file, and then execute the whole sequence by referring to the M-file. Below appears the text for an M-file modelling this dispersal.

%
%     BALLSPORE
%
%     a Matlab M-file which models ballistic dispersal of fungi in 1-D
%

pold = pnew;          % define initial condition, this step --  
                      %     USER MUST DEFINE THE FIRST pnew 
                      %     BEFORE RUNNING THIS PROGRAM

KB=1./(pi*sqrt(c^4/g^2-x.^2)).*(abs(x) <= c^2/g);
                  % define dispersal kernel, taking care that it is zero in 
                  %   the right places

KB=KB/trapz(dx*KB);
                  % this normalizes the probability kernel.  Because of 
                  %   the singularities at the ends this is an important
                  %   step -- otherwise there will be net loss

fKB=dx*fft(KB);   % Take FFT of dispersal kernel, multiplying by dx to save
                  %   that step later

fpnew = fKB.*fft(pold);
                  % FFT of new dispersed population
pnew = real(fftshift(ifft(fpnew)));
                  % new dispersed population

hold on, plot(x,pnew,'r'), hold off
                  % plot results on current figure

When you have this all typed in to the editor, save the file as ballspore.m and click on the command window to activate it. Make sure that the directory you are in in the command window is the same as the directory you save the program ballspore.m in! Now we must set the initial space up, the parameters, and the initial colony:

     g=980;                     % cm/s^2 for gravity constant
     c=50;                      % cm/s speed of ejection
     xlr=15; np=256; dx=2*xlr/np;
     x=-xlr+dx*[0:np-1];        % define independent variable
     pnew=(abs(x) <= 1);

Now we can plot the initial population and call ballspore, which will plot subsequent dispersed populations on the same graph. So, in the command window type:

     plot(x,pnew,'y')
     ballspore

Now, you can see the evolution of the dispersing population by using $\uparrow$ to bring back the ballspore command and then carriage return to execute the command.

Movement at Constant Speed from a Source and Settling

If propagules (density $u$) move at constant speed to the left and right away from a source and settle at some rate $h(t)$, the precipitated individuals (density $v$) satisfy the following system of PDE (from Neubert et. al., 1995):

$$u_t = -c \mbox{ sign} (x) u_x - h(t) u, \quad \quad v_t = h(t) u.$$

If these equations are solved with initial conditions $u(x,0)=\delta(x), v(x,0)=0$ corresponding to a unit release of propagules precisely at the origin, the final distribution of settled individuals ($K_A$) satisfies

$$\frac1{2c} h\left( \frac{\vert x\vert}{c}\right) \exp \left[ -\int_0^{\frac{\vert x\vert}{c}} h(s) ds \right].$$

This solution differs from $K_S$ primarily in that there is no diffusion process operating. Choosing $h(t) = 3 a^3 t^2$ gives an interesting bi-modal solution

$$K_A= \frac{3 a^3 \vert x\vert^2}{2c^3} \exp\left[ -\left(\frac{a \vert x\vert}{c} \right)^3 \right].$$

(see Neubert et. al. for other choices for $h$ and their consequences). Try plotting this kernel for various choices of $c$ and see how it behaves.

EXERCISE 8: Under the assumption that every settled propagule becomes a propagule-producing adult, write an M-file in MATLAB that will allow you to investigate the multi-generational dispersal of an initially-localized population with dispersal probability given by $K_A$. How would you modify this code to reflect mortality or non-germination - say only 50% of dispersed individuals survive to become adults? Could you modify it to have more than one generation producing propagules?

Probability and Dispersal in Two Space Dimensions

Now we have more than enough understanding to move on to two-dimensional dispersal. Mainly there is nothing more complicated than in one dimension, except that it can be computationally much more time consuming and more difficult to visualize. But the computational techniques are identical in MATLAB , except that a few commands (mainly the FFT commands) need to be changed to 2-D versions.

Dispersal via Random Walks in Two Dimensions

As in the case of one dimensional random dispersal, dispersal generated by random walks in two dimensions is modelled by the diffusion equation,

$$u_t = D \nabla^2 u = D (u_{xx} + u_{yy}).$$

The fundamental solution corresponding to initially localized data $u(x,y,0)=\delta(x) \delta(y)$ is given by

$$K_{2D}(x,y,t) = \frac1{4 \pi D t} \exp\left[ -\frac{x^2+y^2}{4 D t} \right].$$

The following MATLAB commands define space and the dispersal kernel K2D:

     xlr=10; ylr=10; np=128; dx=2*xlr/np; dy=2*ylr/np;
     x=linspace(-xlr,xlr-dx,np); 
     y=linspace(-ylr,ylr-dy,np); 
     [X,Y]=meshgrid(x,y);
     t=1; D=1; 
     K2D=1/(4*pi*D*t)*exp(-(X.^2+Y.^2)./(4*D*t));

You may wish to take a second and visualize this using surf or pcolor and shading flat or shading interp. Now, to see how a population initially localized in the rectangle $(-1 \le y \le 1) \times (-3 \le x \le 3)$ disperses, use:

     p0=(abs(X)<=3 & abs(Y)<= 1);
     fp0=fft2(p0);
     fK2D=dx*dy*fft2(K2D);
     fp1=fp0.*fK2D;
     p1=real( fftshift( ifft2( fp1 ) ) );

Some things to observe: the overall format for implementing dispersal is the same as in one dimension, with the exceptions that fft and ifft are replaced by fft2 and ifft2; there is now a multiplication by both dx and dy to implement the convolution. To visualize how the initial field relates to the dispersed field we will use a combination of graphics - one choice would be contours for the initial field and density plots for the dispersed field, as in:

     pcolor(X,Y,p1),shading flat, hold on, contour(X,Y,p0,'m'), hold off

You may want to further use the options axis square and colormap hot. The hot colormap, in particular, is useful in interpreting density plots, since it gives a clear indication of high and low values. To see which colors correspond to which values, you may also use the command colorbar, which will add a color scale with associated scalar values. You might also want to try a combination surface/contour plot, as in

     surf(X,Y,p1+1),shading flat, hold on, contour(X,Y,p0,'m'), hold off

For most applications I find it easiest to use combinations of density and contour plots.

EXERCISE 9: Disperse this population for one more time step, saving p1 and creating a new, twice-dispersed population, p2. Plot all three populations on a single graph using contours in two different colors for p0 and p1 and a hot-colored density plot for p2. Use the combination of commands dx*dy*trapz(trapz(p2)) to estimate the double integral of p2, and then compare that to the previous two dispersal slices to make sure that the population is being conserved.

Modal Dispersal in Two Dimensions

Well, there were many, many different things to do in one dimension, and as you might guess there are infinitely many more in two dimensions. We will discover more of these as we go along, but as a parting shot let's model a combination of advection and diffusion in two dimensions. Suppose a population is performing random walks in a medium which moves in a particular direction with velocities $u$ and $v$ in the $x$ and $y$ directions, respectively. Without belaboring the details, it can be shown that an individual initially localized at the origin disperses with probability

$$K_{2A} = (x,y,t) = \frac1{4 \pi D t} \exp\left[ -\frac{(x-u t)^2+(y-v t)^2}{4 D t} \right].$$

EXERCISE 10: Investigate the advective/diffusive dispersal of a population which is initially localized, according to the probability kernel $K_{2A}$. Pick $D=.25$, $u=2,v=2$, and use the command

     p0=(abs(X+8)<=.5 & abs(Y+8)<=.5);

to specify the initial population. Define the dispersal kernel $K_{2A}$ in MATLAB using the previous grid, and using a sequence of commands on a single line disperse this population over one unit of time. Using a combination of contours and density plots show the population's progress over three time units of dispersal. (Alternatively, you may want to write a single M-file to do this.)

EXERCISE 11: Build an M-file which will allow you to investigate the behavior of a population in an advection/diffusion setting with random winds. That is, over each period of dispersal (`day') the wind speed is constant, but from day to day the direction will vary randomly. Hint: In MATLAB the command for generating uniformly distributed random numbers is rand. Visualize a couple of dispersal steps for a population localized initially uniformly in a disk of radius 1 centered at (0,0). Would it be hard to modify your program to account for random changes in wind speed as well?

Turchin's Model for Prey-Taxis

Another characteristic of moving, living things is that they move for a reason, generally to get somewhere or avoid something. An example would be predators or parasites searching for prey/hosts. Even the simplest of beings is likely to slow down, that is, to take shorter steps in its random searching motions, when it is in the presence of prey items. Let $V(x,y)$ be the density of `victims' in space, and let $P(x,y,t)$ be the density of searchers (predators/parasites). Let $D_2$ be the diffusivity of $P$ in the absence of victims and $D_1 < D_2$ the diffusivity in the presence of victims. Varying motility of the searcher can be modelled (following Turchin, Reeve, Cronin and Wilkens, 1997) by

$$\mu (V) = \frac{D_2 (V + d)}{d + V \frac{D_2}{D_1}} = D_2 - \underbrace{(D_2-D_1) \frac{V}{V+d D_2/D_1}}_{\hat{\mu}},$$

where $d$ is a saturation parameter describing how rapidly motility changes with density of victims. The function $\hat{\mu}$ can be intepreted as the prey-based perturbation to the `normal' motility of predators in search mode. A PDE model for the redistribution of searchers in response to the density of victims is

$$P_t = \nabla^2 \left[ \mu(V) P \right].$$

The idea is that in regions with prey predators slow down (to eat the prey, or at least to hug them and pet them and stroke their fur the wrong way), while in regions without prey they speed up (searching). Thus, when a distribution of predators encounters a distribution of prey the predators tend to `pile up'. In principle this equation is always soluble, but in practice solution requires a different calculation for each choice of $V$. A very approximate solution with initial condition $P_0(x,y)$ is given by:

$$P(x,y,t) = \frac{e^{- \nabla^2 \hat{\mu}(V(x,y)) t}}{4 \pi D... ... (x,y) * \exp\left[ -\frac{x^2+y^2}{4 D_2 t} \right] \right).$$

Remember that $\nabla^2 f = \left(\frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2}\right) f$ and that * denotes the two-dimensional spatial convolution.

This solution becomes increasingly valid as $V$ tends to zero (or as the diffusivities become more alike), and depends critically on $d$ being `large'. The solution is not valid at all for very large or highly variable densities of prey, or when $d$ is small, or when $t$ is taken to be large. Fortunately, solving PDE is not our best and highest goal. As a model for victim-taxis, the approximation above can serve perfectly well! We must separate the elements which give the random character of diffusivity from the elements defining the spatial structure of changing motility in the presence of prey. The approach will be to:

1. convolve with a random motion kernel,
2. multiply by an exponential function ( $e^{-\nabla^2 \hat{\mu} t}$) to give the proper spatial structure, and
3. normalize to ensure that predators have neither been created nor destroyed.

In step 2, an additional crucial operation will be the numerical calculation of $\nabla^2 \hat{\mu}$ using MATLAB's built-in Laplacian calculator (del2). Below appears MATLAB pseudo-code which implements this model for prey-taxis, assuming an initial distribution of prey distributed in a Gaussian oval centered at $x=5=y$,

$$V(x,y) = 100 \exp\left[ -\frac{(x-5)^2}{4} -(y-5)^2 \right],$$

and a population of predators initially distributed uniformly in a circle of radius 3. In particular, notice the `normalization' step which keeps the predator population from growing or shrinking, and the separation between the purely random motion (evaluated with convolution) and the multiplication by the exponential factor, which gives spatial structure to the taxis.

     % 
     %    Testing out Turchin-taxis
     % 
     %    Set up space:
     %
     np=128;
     xl=10; dx=2*xl/np;
     x=linspace(-xl,xl-dx,np);
     [X,Y]=meshgrid(x,x);
     
     %   Define parameters
     
     D2=2; D1=1.5;                     % the motilities
     d=1; t=4;                         % saturation parameter, time
     V=exp(-(X-5).^2./4-(Y-5).^2);     % these are the victims:
   
     %   perturbation to the basic motility (D2):
     muhat=4*del2((D2-D1)*(D2*V)./(d*D1+ D2*V) ,dx,dx);

     %   Set initial population of predators
     
     P=((X.^2+Y.^2)<=4); 
     numP=dx^2*trapz(trapz(P));        % total number predators
     
     %   Define dispersal kernel
     
     K=.25/(pi*D2*t)*exp(-(X.^2+Y.^2)/(4*D2*t));
     fK=dx.^2*fft2(K);  
     
     %   Disperse population according to turchin-taxis
     
     fP=fft2(P);
     Pt=exp(-muhat*t).*real(fftshift(ifft2(fP.*fK)));
     
     %   Normalize the dispersed population
     
     numPt=dx^2*trapz(trapz(Pt));
     Pt=numP/numPt*Pt;
     
     %   Plot the results in comparison with the victims
     
     pcolor(X,Y,Pt), shading interp, colormap hot, axis square
     hold on, contour(X,Y,V,'b'),  hold off

One difficulty with this sort of approach is that it is not `linear' in $t$, the time step. That is, taking ten steps with $t=1$ will be different than taking five steps with $t=2$ or two steps with $t=5$. In general, the approach gets more linear-acting for shorter steps (which is just where one does not wish to use it for an IDE approach). However, it is only a model for chemotaxis, no more or less valid in an absolute sense than the PDE model appearing in the paper of Turchin et. al. (1997). It is important, though, to remember that there is now an extra `parameter', $t$, which must be specified.

EXERCISE 12: With the above MATLAB code for defining a prey-tactic dispersal kernel, determine how a population of predators, initially localized near $x=0=y$, redistributes in response to the population of victims given above. In two different figures plot the density response of the predators as a density field after 4 and 8 time units of dispersal. Indicate the starting point, and plot contours of the victim field to see how the taxis orients the population. Check to see that the number of searching individuals is being preserved. Can you explain these results in terms of decreased motility in the presence of they prey?

Congratulations! You have learned more about probabilistic dispersal and convolutions in MATLAB than most people are learn in a life time!