September 25, 2017, Monday
University of Colorado at Boulder Search A to Z Campus Map University of Colorado at Boulder CU 
Search Links


MBW:Wolf Territoriality, Wolf-Deer Interaction and Survival

From MathBio

Jump to: navigation, search

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

u({\mathbf  {x}},t)={\text{expected density of wolves in pack}},

where {\mathbf  {x}} denotes position and t 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:

J_{u}|_{{{\text{convection}}}}=-\underbrace {\overbrace {c_{u}}^{{{\text{speed of movement}}}}\overbrace {({\mathbf  {x}}-x_{u})}^{{{\text{location of den}}}}}_{{{\text{space-dependent velocity}}}}

A continuous version of c_{u}({\mathbf  {x}}-x_{u}) is:

c_{u}({\mathbf  {x}}-x_{u})=c_{u}\tanh(\beta r){\frac  {{\mathbf  {x}}-x_{u}}{r}},

where r=||{\mathbf  {x}}-x_{u}||. This terms will stop movement of the wolves when they approach their den. c_{u} denotes the maximum speed of the wolf as it moves towards the den, and \beta 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:

J_{u}|_{{{\text{diffusion}}}}=-\overbrace {D(u)}^{{{\text{density-dependent diffusion coefficient}}}}\nabla u.

For our model, we take D(u)=d_{u}u^{n}, where d_{u} is a constant, and n>0. 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:

{\frac  {\partial u}{\partial t}}-\nabla \cdot J_{u}=0\Rightarrow {\frac  {\partial u}{\partial t}}=\nabla \cdot [c_{u}({\mathbf  {x}}-x_{u})u+D_{u}(u)\nabla u]

The specific equation we will use for this model is:

{\frac  {\partial u}{\partial t}}=\nabla \cdot [c_{u}\tanh(\beta r){\frac  {({\mathbf  {x}}-x_{u})}{r}}u+d_{u}u^{n}\nabla u]

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:

J_{n}\cdot {\mathbf  {n}}=0,

where {\mathbf  {n}} is the vector normal to the territory boundary.

For the initial conditions, u({\mathbf  {x}},t=0)=u_{0}({\mathbf  {x}}). In addition, we define the total number of wolves, Q in the pack as:

Q=\int _{\Omega }u({\mathbf  {x}},t)d{\mathbf  {x}}.

where \Omega 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,

{\frac  {\partial }{\partial t}}Q={\frac  {\partial }{\partial t}}\int _{\Omega }u({\mathbf  {x}},t)d{\mathbf  {x}}

0=\int _{\Omega }{\frac  {\partial }{\partial t}}u({\mathbf  {x}},t)d{\mathbf  {x}}

0=-\int _{\Omega }\nabla {\mathbf  {J}}_{u}d{\mathbf  {x}}

0=-\int _{{\partial \Omega }}{\mathbf  {J}}_{u}\cdot {\mathbf  {n}}ds

0=0

We see that the zero-flux boundary conditions agree with our prediction that the Q will not change. We expect that the average population, U_{0} to be: U_{0}={\frac  {1}{A}}\int _{\Omega }u_{0}({\mathbf  {x}})d{\mathbf  {x}}, where A denotes the area of the domain \Omega .

Now let us consider a steady solution to our partial differential equation. That is, let us set {\frac  {\partial }{\partial t}}u({\mathbf  {x}},t)=0. Our equation then becomes:

0=\nabla \cdot [c_{u}({\mathbf  {x}}-x_{u})u+D_{u}(u)\nabla u]

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

c_{u}u\tanh[\beta (x-x_{u})]+d_{u}u^{n}{\frac  {du}{dx}}={\text{constant}}

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

c_{u}u\tanh[\beta (x-x_{u})]+d_{u}{\frac  {du}{dx}}=0

d_{u}{\frac  {du}{dx}}=-c_{u}u\tanh[\beta (x-x_{u})]

{\frac  {du}{u}}=-{\frac  {c_{u}}{d_{u}}}\tanh[\beta (x-x_{u})]dx

\ln u=\ln 1-\ln[\cosh(\beta (x-x_{u}))]^{{{\frac  {c_{u}}{\beta d_{u}}}}}+A,A={\text{constant}}

u={\frac  {e^{A}}{\cosh[\beta (x-x_{u})]^{{{\frac  {c_{u}}{\beta d_{u}}}}}}}

u_{s}(x)={\frac  {B}{\cosh[\beta (x-x_{u})]^{{{\frac  {c_{u}}{\beta d_{u}}}}}}},B=e^{A}.

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

Q=\int _{\Omega }{\frac  {B}{\cosh[\beta (x-x_{u})]^{{{\frac  {c_{u}}{\beta d_{u}}}}}}}

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

c_{u}u\tanh[\beta (x-x_{u})]+d_{u}u^{n}{\frac  {du}{dx}}={\text{constant}}

If we integrate this, we get:

u_{s}(x)={\bigg [}{\frac  {c_{u}n}{d_{u}\beta }}\ln {\bigg (}{\frac  {\cosh(\beta x_{b})}{\cosh(\beta (x-x_{u}))}}{\bigg )}{\bigg ]}^{{{\frac  {1}{n}}}},

where x_{b} denotes the radius of the pack. The total number of wolves is given by

Q=\int _{{x_{u}-x_{b}}}^{{x_{u}+x_{b}}}{\bigg [}{\frac  {c_{u}n}{d_{u}\beta }}\ln {\bigg (}{\frac  {\cosh(\beta x_{b})}{\cosh(\beta (x-x_{u}))}}{\bigg )}{\bigg ]}^{{{\frac  {1}{n}}}}dx.

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: u(x,t)={\text{density of wolf pack 1}}

v(x,t)={\text{density of wolf pack 2}}

p(x,t)={\text{density of RLU markings made by pack 1}}

q(x,t)={\text{density of RLU markings made by pack 2}}

We want our model to reflect the following:

{\text{Rate of change in expected density of pack 1}}

={\text{Rate of change of convective movements of pack 1 wolves towards their den}}

{\text{  }}+{\text{Rate of change due to dispersal of pack 1 wolves}}

{\text{  }}+{\text{Rate of change due to movement of pack 1 wolves away from pack 2 RLUs}}

and

{\text{Rate of change in expected density of pack 2}}

={\text{Rate of change of convective movements of pack 2 wolves towards their den}}

{\text{  }}+{\text{Rate of change due to dispersal of pack 2 wolves}}

{\text{  }}+{\text{Rate of change due to movement of pack 2 wolves away from pack 1 RLUs.}}

We will additionally assume that the two packs are identical in every way, except where their den is located (x_{u}{\text{ and }}x_{v} 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, {\mathbf  {J}}_{c}, which is given by:

{\mathbf  {J}}_{{c_{u}}}=-c_{u}(x-x_{u},q)u

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

c_{u}=C_{u}\tanh(\beta r){\Bigl (}{\frac  {x-x_{u}}{r}}{\Bigr )}{\Bigl (}{\frac  {Aq}{B+q}}{\Bigr )},

where C_{u},\beta ,r 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, {\mathbf  {J}}_{a}, given by:

{\mathbf  {J}}_{{a_{u}}}=a_{u}(q)u\nabla q,

where a_{u}(q) is a monotonically non-decreasing function, much like the one given in the convection flux term just described.

The divergence term, {\mathbf  {J}}_{{d_{u}}} is the same for this model, as it is for the single-pack model, where we take n 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:

{\frac  {\partial p}{\partial t}}=u[l_{p}+m_{p}(q)]-f_{p}p,

where l_{p}{\text{and }}f_{p} are constants describing low-level RLU marking and first-order decay constants, respectively. The m_{p}(q) term describes the the response to foreign RLU markings, and is a bounded and monotonically non-decreasing function, such as the ones described for c_{u} 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:

[{\mathbf  {J}}_{{c_{u}}}+{\mathbf  {J}}_{{d_{u}}}+{\mathbf  {J}}_{{a_{u}}}]\cdot {\mathbf  {n}}=0{\text{ on }}\partial \Omega

and

[{\mathbf  {J}}_{{c_{v}}}+{\mathbf  {J}}_{{d_{v}}}+{\mathbf  {J}}_{{a_{v}}}]\cdot {\mathbf  {n}}=0{\text{ on }}\partial \Omega ,

where \Omega denotes the territory, and {\mathbf  {n}} denotes the vector parallel to the boundary. We also define our initial conditions to be:

u(x,0)=u_{0}(x){\text{, }}v(x,0)=v_{0}(x){\text{, }}p(x,0)=p_{0}(x){\text{, }}q(x,0)=q_{0}(x).

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

\int _{\Omega }u(x,t)dx,

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 A=\int _{\Omega }dx, and the average density of wolves is given by:

{\frac  {1}{U_{0}}}={\frac  {1}{A}}\int _{\Omega }u_{0}(x)dx,{\text{  }}{\frac  {1}{V_{0}}}={\frac  {1}{A}}\int _{\Omega }v_{0}(x).dx

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, L=A^{{{\frac  {1}{m}}}}, 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:

u^{*}={\frac  {u}{U_{0}}},{\text{ }}v^{*}={\frac  {v}{V_{0}}},{\text{ }}p^{*}={\frac  {pf_{p}}{U_{0}l_{p}}},{\text{ }}q^{*}={\frac  {qf_{q}}{V_{0}l_{q}}},{\text{ }}t^{*}=tf_{p},x^{*}={\frac  {x}{L}},{\text{ }}c_{u}^{*}={\frac  {c_{u}}{Lf_{p}}},

c_{v}^{*}={\frac  {c_{v}}{Lf_{q}}},{\text{ }}d_{u}^{*}={\frac  {d_{u}}{L^{2}f_{p}}},{\text{ }}d_{v}^{*}={\frac  {d_{v}}{L^{2}f_{q}}},{\text{ }}a_{u}^{*}={\frac  {a_{u}V_{0}l_{q}}{L^{2}f_{p}^{2}}},{\text{ }}a_{v}^{*}={\frac  {a_{v}U_{0}l_{p}}{L^{2}f_{q}^{2}}}

m_{p}^{*}={\frac  {m_{p}}{l_{p}}},{\text{ }}m_{q}^{*}={\frac  {m_{q}}{l_{q}}},{\text{ }}\phi ={\frac  {f_{p}}{f_{q}}}

The non dimensionalized system now becomes:

{\frac  {\partial u}{\partial t}}+\nabla \cdot [{\mathbf  {J}}_{{c_{u}}}+{\mathbf  {J}}_{{d_{u}}}+{\mathbf  {J}}_{{a_{u}}}]=0

{\frac  {\partial v}{\partial t}}+\nabla \cdot [{\mathbf  {J}}_{{c_{v}}}+{\mathbf  {J}}_{{d_{v}}}+{\mathbf  {J}}_{{a_{v}}}]=0

{\frac  {\partial p}{\partial t}}=u[1+m_{p}(q)]-p

{\frac  {\partial q}{\partial t}}=v[1+m_{q}(p)]-\phi q

where

{\mathbf  {J}}_{{c_{u}}}=-uc_{u}(x-x_{u},q),{\text{ }}{\mathbf  {J}}_{{d_{u}}}=-d_{u}(u)\nabla u,{\text{ }}{\mathbf  {J}}_{{a_{u}}}=a_{u}(q)u\nabla q

{\mathbf  {J}}_{{c_{v}}}=-vc_{v}(x-x_{v},p),{\text{ }}{\mathbf  {J}}_{{d_{v}}}=-d_{v}(v)\nabla v,{\text{ }}{\mathbf  {J}}_{{a_{v}}}=a_{v}(q)v\nabla p.

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

u_{0}^{*}={\frac  {u_{0}}{U_{0}}},{\text{ }}v_{0}^{*}={\frac  {v_{0}}{V_{0}}},{\text{ }}p_{0}^{*}={\frac  {p_{0}f_{p}}{U_{0}l_{p}}},{\text{ }}q_{0}^{*}={\frac  {q_{0}f_{q}}{V_{0}l_{q}}}.

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 \Omega 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:

{\frac  {\partial u}{\partial t}}=\nabla \cdot [c_{u}(x-x_{u},q)u+d_{u}\nabla u]

{\frac  {\partial v}{\partial t}}=\nabla \cdot [c_{v}(x-x_{v},p)v+d_{v}\nabla v]

{\frac  {\partial p}{\partial t}}=u[1+m_{p}(q)]-p

{\frac  {\partial q}{\partial t}}=v[1+m_{q}(p)]-q.

For c_{u}{\text{ and }}c_{v}, we use the formula:

c_{u}=C_{u}\tanh(\beta r){\Bigl (}{\frac  {x-x_{u}}{r}}{\Bigr )}{\Bigl (}{\frac  {Aq}{B+q}}{\Bigr )}.

We use a similar formula for c_{v}. For m_{p}(q), we use the formula:

m_{p}(q)={\frac  {Aq}{B+q}}

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:

{\frac  {\partial u}{\partial t}}=\nabla \cdot [c_{u}(x-x_{u},q)u+d_{u}\nabla u-a_{u}(q)u\nabla q]

{\frac  {\partial v}{\partial t}}=\nabla \cdot [c_{v}(x-x_{v},p)v+d_{v}\nabla v-a_{v}(p)v\nabla p]

{\frac  {\partial p}{\partial t}}=u[1+m_{p}(q)]-p

{\frac  {\partial q}{\partial t}}=v[1+m_{q}(p)]-\phi q,

where m_{p}(q),{\text{ }}m_{q}(p),{\text{ }}a_{u}(q),{\text{ }}a_{v}(p),{\text{ }}c_{u}(x-x_{u},q),{\text{ and }}c_{v}(x-x_{v},p) 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')