This is a summary of Chapter 14, Sections 1-3, in *Mathematical Biology II: Spatial Models and Biomedical Applications* by J.D. Murray. This summary was done by Takako Hirokawa.

## Introduction

How predators and their prey coexist in the same region, especially when predators are territorial creatures, is something that is not very well understood. Additionally, how are these territories determined and how are they maintained? Wolves are territorial creatures that have been well studied over the last several decades and therefore provide a wealth of knowledge that can help us understand how territorial animals behave. Most of the current data on wolves comes from radio-tagged wolves from northeastern Minnesota, where a large portion of information comes from Isle Royale National Park, an isolated wolf-moose ecosystem.

Wolves are social creatures that live in very strict social structures called packs. These packs consist of 3-15 animals and are lead by a pair of so-called alpha wolves. They form territories that are related to the size of their pack. These territories rarely overlap, and are sometimes maintained over several years. In addition, between the territories, there are buffer zones of about 2 km. These buffer zones are like 'no man's lands'. This is further supported by the behavior of the wolves, who seem to stay strictly within pack territory boundaries. Territories are maintained through scent-marking and howling. Raised leg urination (RLU) is the most common form of scent-marking. These markings increase around the territory boundaries. It has been observed that deer and moose, the wolf's main choice of prey, stay generally in the buffer zones. Territorial zones result from the steady state solution to the model equations.

## Models for Wolf Pack Territory Formation: Single Pack

In this section, we will focus simply on a single pack and how the population is distributed within the territory. With our model, we would like to be able to answer the following questions:

- Can we determine how pack territories form as well as determine their size and why the territories are stable for many years?

- Can we explain why there are more conflicts between wolves, buffer trespassing, and territory changes during the winter months?

- Is there a way to predict wolf behavior during times of low prey populations?

- Can we quantify these predictions based on behavioral parameters?

- What role do the buffer zones play in wolf-deer interactions? Do they stabilize both populations because the deer have a refuge in these buffer zones and are therefore prevented from being extinct?

- As both wolves and deer are affected by seasonal changes in each of their populations, how do migrations affect wolf-deer interactions?

First, we will focus on a single pack territory and how the population is distributed within the territory. Let

where denotes position and denotes time.

Our model for the single pack will be that the rate of change in wolf density is equal to the rate of change due to movement of the wolves towards the den, the social center for the wolves, in addition to the rate of change due to dispersal of wolves in search of food.

To determine the rate due to movement towards the den, we will assume that the movements are more centered around the den during the summer months due to puppy-rearing, and that movements toward the den from anywhere within the territory are generally in a straight line. Interestingly, this last point is in fact exhibited by the wolves who are well aware of where they located within their territory.

In general, this rate term can be modeled as a convection term:

A continuous version of is:

where . This terms will stop movement of the wolves when they approach their den. denotes the maximum speed of the wolf as it moves towards the den, and denotes the change in the rate of convective movement.

The second term, the rate of change in the density due to food scavenging, we make two assumptions. The first is that food is plentiful and homogeneously distributed throughout the territory. This assumption leads to our second; the wolves do not have a biased direction of movement because the food is uniformly distributed. This second term can be modeled generally as a diffusion flux term:

For our model, we take , where is a constant, and . With a diffusion coefficient such as the one we have chosen, we can model a population that tends stay relatively close to the den, where the wolves are most familiar and the population has the highest density.

Therefore, if we combine the terms we have just derived, we get our general model equation:

The specific equation we will use for this model is:

Now, we can determine the appropriate initial and boundary conditions. For the boundary conditions, we can assume that there are no wolves leaving or entering the territory, as wolves generally stay within territory boundaries. In other words, we will assume zero-flux boundaries:

where is the vector normal to the territory boundary.

For the initial conditions, . In addition, we define the total number of wolves, in the pack as:

where denotes the domain, or the territory.

We expect that the total number of wolves will not change at any given point in time, neglecting birth and death rates within the pack. So,

We see that the zero-flux boundary conditions agree with our prediction that the will not change. We expect that the average population, to be:
where denotes the area of the domain

Now let us consider a steady solution to our partial differential equation. That is, let us set . Our equation then becomes:

We now consider the 1D case, and substitute in our convection term for the continuous version:

For the linear diffusion case, we assume that so the equation above simplifies to:

This last equation is the steady state solution to the differential equation. The total number of wolves in this case is given by:

For the steady state with a nonlinear diffusion term, we start with the differential equation:

If we integrate this, we get:

where denotes the radius of the pack. The total number of wolves is given by

Linear Diffusion term (n=0), x-axis represents distance, y-axis represents u

Nonlinear Diffusion term (n>0), x-axis represents distance, y-axis represents u

We see that when there n=0, we do not get a defined territory boundary, as the number of wolves at the edges of the territory never quite equals zero. This does not agree with the finding that the wolves have very well-defined territories. It is clear that the steady state for the n>0 diffusion term does in fact give definite territory boundaries.

This figure shows the population evolving over time from the initial condition where there are wolves spread evenly over all of the territory and slowly evolves into the steady state.

## Multi-Pack Territory Formation

Before one can start to model several packs, it is important to assume how the wolf packs will interact with each other. Wolf packs respond to a foreign pack's raised leg urination (RLU), but how they react to a foreign pack's RLU is poorly understood. This can be modeled in two different ways. In the first model, the wolves respond to the presence of RLUs by an increased rate of movement back to the den while producing their own RLU markings. In the second, the wolves always want to move down the gradient of density of foreign RLU markings.

In this section, we consider two packs in one dimension. The notation we will use will be the following:

We want our model to reflect the following:

and

We will additionally assume that the two packs are identical in every way, except where their den is located ( for packs 1 and 2, respectively), so that from now on, we can discuss only pack 1 and assume the same is true for pack 2 unless states otherwise. The first term in each of these equations can be modeled by a convection term, , which is given by:

This term assumes that the wolves will respond to foreign RLU markings, , (this would be for pack 2) by increasing both their speed towards the den. In addition, we desire that such that we elicit a response from the wolves. The function should be a bounded, monotonically non-decreasing function in so that the elicited response is greater for a larger . A function such as:

where are as defined above, and A, B are constants.

The other way in which we can model pack 1's response to the RLUs is to include another gradient term, , given by:

where is a monotonically non-decreasing function, much like the one given in the convection flux term just described.

The divergence term, is the same for this model, as it is for the single-pack model, where we take to equal 1.

Now, we must better define the RLU functions. Most of our assumptions come from what has been observed in the pack behavior in northern Minnesota. It has been found that there is always at least some low level of RLU markings, and that the rate of these markings increase considerably when there are foreign RLU markings nearby. These RLU markings also decay over time. Though the rate of decay is highly dependent on environmental conditions, we will assume first-order decay. We can now construct an equation to describe the RLUs:

where are constants describing low-level RLU marking and first-order decay constants, respectively. The term describes the the response to foreign RLU markings, and is a bounded and monotonically non-decreasing function, such as the ones described for above.

Now that we have determined the equations that describe the population density, we must now determine boundary and initial conditions. Just as we had for the single-pack case, we assume that there is neither emigration from nor immigration into the territory. This would translate to no-flux boundary conditions:

and

where denotes the territory, and denotes the vector parallel to the boundary. We also define our initial conditions to be:

Furthermore, we say that the total number of wolves from pack 1 is given by:

which we can find using a similar argument to that we used for the single-pack model, that this agrees with our boundary conditions.

Additionally, the area of the territory is given by , and the average density of wolves is given by:

As this is a very complicated system, we now take the time to nondimensionalise the system, including the initial and boundary conditions. We start by defining a new length, where m denotes the number of dimensions of the territory (which therefore equals 1 or 2). We can nondimensionalise the system in the following way:

The non dimensionalized system now becomes:

where

Note that we have dropped the asterisks for notational simplicity. The boundary conditions are now nondimensionalised to become

From now, we will drop the asterisks, as we did with the model equations. Our initial conditions remain unchanged while, the area of our territory has become 1.

We now state the two ways in which the wolves respond to foreign RLU markings more explicitly. We first consider the assumption that coming into contact with a foreign RLU marking causes the wolf to increase their speed back to their den and to increase the rate of RLU production. Therefore, we will specifically use the following set of equations:

For , we use the formula:

We use a similar formula for . For , we use the formula:

Wolf Pack 1, represented in the equations as u.

Wolf Pack 2, represented in the equations as v.

We see that from the two plots of the wolf packs, that the two packs seems to concentrate themselves in two different areas of the domain, hence creating territories.

For the second type of RLU marking, we assume that the wolves respond to a gradient of foreign RLU markings, so we modify our set of modeling equations (and nondimensionalise) to get:

where are given above.

## References

*Mathematical Biology II: Spatial Models and Biomedical Applications* by J.D. Murray

"On Wolf Territoriality and Deer Survival" by K.A. Jane White, M.A. Lewis, and J.D. Murray, found in *Modeling Spatiotemporal *

"Wolf-Deer Interactions: A Mathematical Model" by K.A.J. White, J.D. Murray, and M.A. Lewis, found in *Proceedings: Biological Sciences*, Vol. 263, No. 1368

"Modelling territoriality and wolf-deer interactions" by M.A. Lewis, and J.D. Murray, found in *Letters to Nature*, Vol. 366,

## MATLAB Code

This m-file generates the contour plot that shows a single wolf-pack evolving over time.

% Toy Model: Evolution of the Wolf Population -- 1 D
% Code by Takako Hirokawa
% 3/31/2011
clear all
close all
digits(50)
% Variables:
% u = spatial distribution - vector, depends on x and t
% x = spatial vector
% us = theoretical steady state solution
% Q = total number of wolves
%%%%%%%%%%%%%%%%%%%% PARAMETERS %%%%%%%%%%%%%%%%%%%%%%
N = 100; % Number of points along x
Nt = 50000; % number of time iterations
beta = .5; % change in rate of convective movement
n = 1; % part of the diffusion coeff
cu = .5; % max speed of wolf when moving towards the den
du = .03; % diffusion coefficient
a = 0; % left end
b = 10; % right end
deltat = .001; % time interval
deltax = (b-a)/N; % space interval
xu = a+(b-a)/2; % location of den -- halfway between a & b
divisor = 1000;
Tt= Nt/divisor;
k=1;
% Construct vectors
t = zeros(1,Nt); % time vector
plott = zeros(1,Tt);
for i = 1:Nt
t(i) = i*deltat;
if mod(i,divisor)==0
plott(k) = t(i);
k = k+1;
end
end
x = zeros(N,1); % spatial vector
for i = 1:N
x(i) = deltax*i+a;
end
uinitial = zeros(N,1); % initial conditions for each x
for i = 1:N
uinitial(i) = 1;
end
u = zeros(N,Nt); % population density array
u(:,1) = uinitial(:,1); % set the initial conditions to the first column
plotmat = zeros(N,Tt);
plotmat(:,1) = uinitial(:,1);
k = 2;
% Use finite difference equation for the PDE
for m = 1:Nt-1
u(1,m+1) = deltat*(du*(u(1,m)^n)*(u(2,m)-2*u(1,m))/(deltax^2)+cu*u(1,m))+u(1,m);
for j = 2:N-1
if x(j) == xu
u(j,m+1) = deltat*(cu*beta*sech(beta*abs(x(j)-xu))^2*((u(j+1,m)-u(j-1,m))/(2*deltax))+du*(u(j,m)^n) ...
*(u(j-1,m)-2*u(j,m)+u(j+1,m))/(deltax^2)+cu*u(j,m))+u(j,m);
else
u(j,m+1) = deltat*(cu*beta*sech(beta*abs(x(j)-xu))^2*((u(j+1,m)-u(j-1,m))/(2*deltax))+du*(u(j,m)^n) ...
*(u(j-1,m)-2*u(j,m)+u(j+1,m))/(deltax^2)+cu*u(j,m))+u(j,m);
end
end
u(N,m+1) = deltat*(du*(u(N,m)^n)*(u(N-1,m)-2*u(N,m))/(deltax^2)+cu*u(N,m))+u(N,m);
if mod(m,divisor)==0
plotmat(:,k) = u(:,m);
k=k+1;
end
end
surf(plott,x,plotmat)
xlabel('t')
ylabel('x')

The following code gives two plots of two wolf packs, in one dimension.

% Toy Model: Evolution of Multiple Wolf Populations -- 1 D
% Code by Takako Hirokawa
% Modified: 4/30/2011
clear all
close all
digits(50)
%%%%%%%%%%%%%%%%%%%% PARAMETERS/CONSTANTS %%%%%%%%%%%%%%%%%%%%%%
N = 200; % Number of points along x
Nt = 50000; % number of time iterations
beta = .5; % change in rate of convective movement
n = 1; % part of the diffusion coeff
cu = 3; % max speed of wolf when moving towards the den
du = .03; % diffusion coefficient
cv = 3; % max speed of wolf when moving towards the den
dv = .03; % diffusion coefficient
a = 0; % left end
b = 1; % right end
deltat = .0001; % time interval
deltax = (b-a)/N; % space interval
xu = 0; % location of den for pack 1
xv = 1; % location of den for pack 2
divisor = 1000; % divisor
A = 1; % constant in J_cu, J_cv, J_au, J_av, m_p, m_q
B = 10; % constant in J_cu, J_cv, J_au, J_av, m_p, m_q
Tt= Nt/divisor;
k=1;
% Construct vectors
t = zeros(1,Nt); % time vector
plott = zeros(1,Tt);
for i = 1:Nt
t(i) = i*deltat;
if mod(i,divisor)==0
plott(k) = t(i);
k = k+1;
end
end
x = zeros(N,1); % spatial vector
for i = 1:N
x(i) = deltax*i+a;
end
uinitial = ones(N,1); % initial conditions for each x
vinitial = ones(N,1);
pinitial = ones(N,1);
qinitial = ones(N,1);
for i = 1:N
uinitial(i) = 1/N;
vinitial(i) = 1/N;
pinitial(i) = 1/N;
qinitial(i) = 1/N;
end
% Construct matrices
u = zeros(N,Nt); % population density array for pack 1
v = zeros(N,Nt); % population density array for pack 2
p = zeros(N,Nt); % RLU density for pack 1
q = zeros(N,Nt); % RLU density for pack 2
% Set initial conditions
u(:,1) = uinitial; % set the initial conditions to the first column for pack 1
v(:,1) = vinitial; % ditto for pack 2
p(:,1) = pinitial; % set the initial conditions for RLU for pack 1
q(:,1) = qinitial; % ditto for pack 2
plotu = zeros(N,Tt);
plotu(:,1) = uinitial;
plotv = zeros(N,Tt);
plotv(:,1) = vinitial;
k = 2;
% Use finite difference equation for the PDE
for m = 1:Nt-1
u(1,m+1) = deltat*(cu*beta*(sech(beta*abs(x(1)-xu)))^2*sign(x(1)-xu)*(A*q(1,m)/(B+q(1,m)))*u(1,m) - ...
du*u(1,m)*((u(2,m)-u(1,m))/(deltax^2)))+u(1,m);
p(1,m+1) = u(1,m)*deltat*(1+A*q(1,m)/(B+q(1,m)))+p(1,m)*(1-deltat);
v(1,m+1) = deltat*(cv*beta*(sech(beta*abs(x(1)-xv)))^2*sign(x(1)-xv)*(A*p(1,m)/(B+p(1,m)))*v(1,m) - ...
dv*v(1,m)*((v(2,m)-v(1,m))/(deltax^2)))+v(1,m);
q(1,m+1) = v(1,m)*deltat*(1+A*p(1,m)/(B+p(1,m)))+q(1,m)*(1-deltat);
for j = 2:N-1
u(j,m+1)=deltat*(cu*beta*(sech(beta*abs(x(j)-xu)))^2*sign(x(j)-xu)*(A*q(j,m)/(B+q(j,m)))*u(j,m)+ ...
cu*tanh(beta*abs(x(j)-xu))*(B/((B+q(j,m))^2))*((q(j+1,m)-q(j-1,m))/(2*deltax))*u(j,m)+ ...
cu*tanh(beta*abs(x(j)-xu))*(A*q(j,m)/(B+q(j,m)))*((u(j+1,m)-u(j-1,m))/(2*deltax))+ ...
du*((u(j+1,m)-2*u(j,m)+u(j-1,m))/(deltax^2)))+u(j,m);
p(j,m+1) = u(j,m)*deltat*(1+A*q(j,m)/(B+q(j,m)))+p(j,m)*(1-deltat);
v(j,m+1)=deltat*(cv*beta*(sech(beta*abs(x(j)-xv)))^2*sign(x(j)-xv)*(A*p(j,m)/(B+p(j,m)))*v(j,m)+ ...
cv*tanh(beta*abs(x(j)-xv))*(B/((B+p(j,m))^2))*((p(j+1,m)-p(j-1,m))/(2*deltax))*v(j,m)+ ...
cv*tanh(beta*abs(x(j)-xv))*(A*p(j,m)/(B+p(j,m)))*((v(j+1,m)-v(j-1,m))/(2*deltax))+ ...
dv*((v(j+1,m)-2*v(j,m)+v(j-1,m))/(deltax^2)))+v(j,m);
q(j,m+1) = v(j,m)*deltat*(1+A*p(j,m)/(B+p(j,m)))+q(j,m)*(1-deltat);
end
u(N,m+1) = deltat*(cu*beta*(sech(beta*abs(x(N)-xu)))^2*sign(x(N)-xu)*(A*q(N,m)/(B+q(N,m)))*u(N,m) - ...
du*u(j,m)*((-u(N,m)+u(N-1,m))/(deltax^2)))+u(N,m);
p(N,m+1) = u(N,m)*deltat*(1+A*q(N,m)/(B+q(N,m)))+p(N,m)*(1-deltat);
v(N,m+1) = deltat*(cv*beta*(sech(beta*abs(x(N)-xv)))^2*sign(x(N)-xv)*(A*p(N,m)/(B+p(N,m)))*v(N,m) - ...
dv*v(j,m)*((-v(N,m)+v(N-1,m))/(deltax^2)))+v(N,m);
q(N,m+1) = v(N,m)*deltat*(1+A*p(N,m)/(B+p(N,m)))+q(N,m)*(1-deltat);
if mod(m,divisor)==0
plotu(:,k) = u(:,m);
plotv(:,k) = v(:,m);
k=k+1;
end
end
mesh(plott,x,plotu)
xlabel('t')
ylabel('x')
title('Wolf Pack 1')
figure
mesh(plott,x,plotv)
xlabel('t')
ylabel('x')
title('Wolf Pack 2')