%\section{Modelling Approach}
%
%We use the discussion in the previous section to justify
%our quantitative modelling approach. We propose using
%two types of models, the first of which:
%\begin{enumerate}
%\item has a time horizon of several months,
%\item is a spatio-temporal model which can be used to analyze the
%formation of territorial structures and resulting spatial heterogeneities,
%\item stresses the behavioral aspects of the
%wolf-wolf and wolf-deer interactions (as opposed to
%the ecological aspects, which typically
%take place over years),
%\item employs a continuum approach to describing expected densities of
%wolves and deer,
%\item takes the season (summer or winter) into account,
%\item does not include birth (because this happens once a year), but
%includes mortality (because this is a continuous process).
%\end{enumerate}
%Such continuum models can be shown to be the result using behavioral `rules'
%in a `random walk' model for individuals (see also Section 5.1).
%The second model
%\begin{enumerate}
%\item has a time horizon of many years,
%\item incorporates seasonal events, such as the spring-time wolf
%and deer births,
%and deer migrations in the spring and fall.
%\item can be used to analyze the degree to which spatial refuges and migration
%help to stabilize wolf-deer interactions
%\item stresses the ecological aspects of the wolf-deer interactions,
%\item employs a discrete, difference equation approach to describing the
%\item includes juvenile and adult classes for wolves and deer (see Hastings, 1984).
%populations of wolves and deer.
%\end{enumerate}
%To distinguish between these two models, we label the first as the {\em behavioral}
%model and the second as the {\em ecological} mode because these features are
%emphasized in the respective models. In actual fact, both models have behavioral
%and ecological components.
\section{Models for Territorial Pattern Formation in Wolves}
In this section we derive mathematical equations
to describe the behavioral interactions between two adjacent
wolf packs. These interactions play a crucial role in developing and
maintaining territory and in regulating the deer populations.
We give the model equations in words prior to their mathematical formulation.
Due to the small number of wolves in a particular pack,
it is likely that, at any given time, wolves may not be found
on large regions of their home territory. Thus, it makes the most
sense to employ a probabilistic approach and
consider the {\em expected} densities of wolves at
point $\ul{x}$ and time $t$, as direct field observations typically will
not yield the exact densities.
Thus, the pertinent state variables are:
\begin{eqnarray}
u(\ul{x},t)&=&\mbox{expected density of wolves from pack number 1}\nonumber\\
v(\ul{x},t)&=&\mbox{expected density of wolves from pack number 2}\nonumber\\
p(\ul{x},t)&=&\mbox{expected density of RLUs from wolf pack number 1}\nonumber\\
q(\ul{x},t)&=&\mbox{expected density of RLUs from wolf pack number 2}\nonumber
\end{eqnarray}
\subsection{Model for a Single Pack}
In the simplest possible case, with only
one wolf pack, we anticipate that two types of movements will dominate:
\begin{enumerate}
\item dispersal as the wolves search for food and other resources,
and
\item movement back towards the den as the wolves return to the social
organizing center (summer only).
\end{enumerate}
Thus the word equation is:
\begin{eqnarray}
\lefteqn{\mbox{Rate of change in expected wolf density}}\nonumber\\
&=&
\mbox{Rate of change due to movement of wolves towards the den (summer only)}\nonumber\\
&+&
\mbox{Rate of change due to dispersal of wolves away from high density} \nonumber\\
& &
\mbox{regions in search of food and other resources.} \nonumber
\end{eqnarray}
We model this situation by a nonlinear density-dependent diffusion-convection
equation of the following form:
\begin{equation}
\tpa{u}=\nabla\cdot\left[c_u
\frac{(\ul{x}-\ul{x}_u)}{|\ul{x}-\ul{x}_u|}
u\right]
+\nabla\cdot\left[d_u(u)\nabla u\right]. \label{eq:w1}
\end{equation}
The convection speed towards the den is scaled by $c_u$;
the magnitude of the diffusion away
from the den is given by the nonlinear function $d_u(u)$;
and the location of the den is $\ul{x}_u$.
A related equation was proposed by Okubo
(1986) to describe the swarming of insects ({\em c.f.} Murray, 1989).
\subsection{Analytical Solutions for the Single-Pack Model}
The simple single-wolf pack model (\ref{eq:w1}) is advantageous because
analytical solutions are known for a number of special cases.
It is, for example, natural to expect that when opposing diffusive and
convective fluxes balance, there may be a spatially heterogeneous steady state
solution. On an infinite solution
domain, such a one-dimensional steady state solution can be calculated
exactly for a special case where
\begin{equation}
d_u(u)=d_0u^\alpha. \label{eq:nlfn}
\end{equation}
Integrating (\ref{eq:w1})--(\ref{eq:nlfn}) twice using the
conditions $u\rightarrow 0$,
$\partial u/\partial x \rightarrow 0$ as $|x|\rightarrow\infty$ yields
the steady state solution
\begin{equation}
u(x)=\left\{
\begin{array}{ll}
\left(\frac{c_u}{d_0}\right)^{1/\alpha}\left[\frac{U_0(1+\alpha)}{2}-\alpha |x-x_u|\right]^{1/\alpha},
&|x-x_u|\leq \frac{U_0(1+\alpha)}{2\alpha}\\[3ex]
0,&|x-x_u|>\frac{U_0(1+\alpha)}{2\alpha}
\end{array}
\right.
\label{eq:w1_sol}
\end{equation}
($\alpha>0$) where $U_0$ is the total number of wolves
\begin{displaymath}
U_0=\int_{-\infty}^{\infty} u(x)\, dx
\end{displaymath}
(Figure 1) ({\em c.f.} Murray, 1989).
When $\alpha=0$ the solution is
\begin{equation}
u(x)=\left\{
\begin{array}{ll}
\left(\frac{U_0c_u}{d_0}\right)^{1/2}- \frac{c_u |x-x_u|}{d_0},
&|x-x_u|\leq \left(\frac{U_0d_0}{c_u}\right)^{1/2}\\[3ex]
0,&|x-x_u|>\left(\frac{U_0d_0}{c_u}\right)^{1/2}.
\end{array}
\right.
\label{eq:w1_sol_a0}
\end{equation}
Based on our discussion elsewhere, $c_u=0$ during the winter when the
den or rendezvous site is no longer the focus of social activity.
When $\alpha=0$ (\ref{eq:w1})--(\ref{eq:nlfn}) becomes a simple diffusion
equation whose solution can be calculated exactly in terms of the Green's
function for the diffusion operator. In one spatial dimension
the solution is
\begin{equation}
u(x,t)=\int_{-\infty}^{\infty}\frac{u(y,0)}{2(\pi d_0 t)^{1/2}}
\exp\left(\frac{-(x-x_u-y)^2}{4d_0t}\right)\,dy.
\label{eq:dif_sol}
\end{equation}
Thus, for example, if all the wolves are initially at the den, then
\begin{equation}
u(\ul{x},0)=U_0\delta(\ul{x}-\ul{x}_u)
\label{eq:delta_ic}
\end{equation}
and (\ref{eq:dif_sol}) describes a radially
expanding, normally distributed population.
When $c_u=0$ and $\alpha>0$ (\ref{eq:w1})--(\ref{eq:nlfn})
becomes the porous media equation.
Weak solutions for the radially symmetric case
with the initial condition (\ref{eq:delta_ic})
can be calculated exactly by using a similarity transformation
(see, for example, Murray (1989) or Grindrod (1991)).
In the one-dimensional case, the solution,
\begin{equation}
u(x,t)=\left\{
\begin{array}{ll}
\frac{1}{\lambda(t)}\left[1-\left(\frac{x-x_u}{r_0\lambda(t)}\right)^2
\right]^{1/\alpha},&|x-x_u|\leq r_0\lambda(t)\\[3ex]
0,&|x-x_u|>r_0\lambda(t),
\end{array}
\right.
\label{eq:pm_1d}
\end{equation}
where
\begin{displaymath}
\lambda(t)=\left(\frac{t}{t_0}\right)^{1/(2+\alpha)},\;\;\;
r_0=\frac{U_0\Gamma(\frac 1 \alpha + \frac 3 2)}
{\pi^{1/2}\Gamma(\frac 1 \alpha + 1)},\;\;\;
t_0=\frac{r_0^2\alpha}{2d_0(\alpha+2)},
\end{displaymath}
($\Gamma$ is the gamma function)
describes a radially expanding wave (Figure 1) (Murray, 1989).
However, the simple situation, described above,
considers neither the population dynamics for the prey species
(white-tailed deer) nor the competition between adjacent packs.
Thus the advantage of being able to calculate exact solutions analytically
is offset by the fact that important interactions are ignored.
We now describe a model which includes these dynamics.
\subsection{Wolf-Wolf Territorial Model}
Inter-pack competition and subsequent territorial behavior
is modelled by
\begin{enumerate}
\item changing equation (\ref{eq:w1}) to include
responses both familiar RLUs and to foreign RLUs from
other packs, and
\item by including equations for the RLU densities
themselves.
\end{enumerate}
We assume that when members of a pack encounter RLUs from
an adjacent pack, they tend to move away from these foreign RLUs and
back towards the den while also increasing their rate of RLU marking.
When pack members encounter RLUs from their own pack, they move towards
these familiar RLUs, possibly following a trail of RLU marks.
Although mortal strife may occur when adjacent packs interact,
for the purpose of modelling the populations
we assume that such fatal interactions are very rare and that the number of
wolves remains effectively constant over the time horizon for the model.
We consider the case where there are two competing
packs, with probability densities $u$ and $v$. Thus
we have:
\begin{eqnarray*}
\lefteqn{\mbox{Rate of change in expected density of wolves (pack 1)}}\nonumber\\
&=&
\mbox{Rate of change due to movement of pack 1 wolves towards their den}\nonumber\\
&+&
\mbox{Rate of change due to dispersal of pack 1 wolves}\nonumber\\
&+&
\mbox{Rate of change due to movement of pack 1 wolves away from the RLUs}\nonumber\\
& &\mbox{made by pack 2}\nonumber\\
&+&
\mbox{Rate of change due to movement of pack 1 wolves towards the RLUs}\nonumber\\
& &\mbox{made by pack 1.}\nonumber\\
\end{eqnarray*}
\begin{equation}
\tpa{u}=\nabla\cdot\left\{\ul{J}_{c_u}+\ul{J}_{d_u}+\ul{J}_{a_u}+\ul{J}_{b_u}
\right\}
\label{eq:w1a}
\end{equation}
where the fluxes are given by
\begin{displaymath}
\begin{array}{rclcc}
\ul{J}_{c_u}&=&u c_u(q) \frac{(\ul{x}-\ul{x}_u)}{|\ul{x}-\ul{x}_u|},&
c_u(0)\geq 0,&\zd{c_u}{q}\geq 0\\[2ex]
\ul{J}_{d_u}&=&d_u(u) \nabla u,&
d_u(0)\geq 0,&\zd{d_u}{u}\geq 0\\[2ex]
\ul{J}_{a_u}&=&ua_u(q) \nabla q,&
a_u(0)\geq 0,&\zd{a_u}{q}\geq 0\\[2ex]
\ul{J}_{b_u}&=&-ub_u(p) \nabla p,&
b_u(0)\geq 0.& \\[2ex]
\end{array}
\end{displaymath}
Here the rate of movement back to the den (which is situated at $\ul{x}_u$)
may depend on the expected levels of competing RLUs ($q$). In the presence
of competing RLUs, the wolves may tend to retreat towards the den and thus
$d c_u/dq\geq 0$. Also, in the spring and summer, when pack movements
are focused upon nourishing the pups at the den or rendezvous site, there may
be an additional constant rate of movement towards $\ul{x}_u$ ($c_u(0)>0$).
The diffusive flux ($d_u(u)$) is similar to that described in the earlier
simpler model (\ref{eq:w1}).
Formulations for the fluxes away from competing
RLUs ($\ul{J}_{a_u}$) and towards familiar RLUs ($\ul{J}_{b_u}$)
assume that movement is in response to the gradient in expected RLU density
(either $\nabla q$ or $\nabla p$). The nonlinear dependence of this
movement upon expected RLU densities (either $q$ or $p$) is described
by $a_u(q)$ and $b_v(p)$. Note that the different signs for $\ul{J}_{a_u}$
and $\ul{J}_{b_u}$ reflect that $\ul{J}_{a_u}$ describes movement {\em away}
from competing RLUs and $\ul{J}_{b_u}$ describes movement {\em towards} familiar
RLUs. Similar formulations for taxes (movements in response to
gradients) have been used to describe the movement of bacteria up
chemical gradients (see, for example, Keller and Segel, date), the
movement of predators up gradients in prey density (Kareiva and Odell, date),
and the movement of herbivores up gradients in plant density (Gueron and Liron,
date) or in plant quality (Lewis, submitted).
The equation for movement of the second wolf pack mirrors that
describing movement of wolf pack 1:
\begin{eqnarray*}
\lefteqn{\mbox{Rate of change in expected density of wolves (pack 2)}}\nonumber\\
&=&
\mbox{Rate of change due to movement of pack 2 wolves towards their den}\nonumber\\
&+&
\mbox{Rate of change due to dispersal of pack 2 wolves}\nonumber\\
&+&
\mbox{Rate of change due to movement of pack 2 wolves away from the RLUs}\nonumber\\
& &\mbox{made by pack 1}\nonumber\\
&+&
\mbox{Rate of change due to movement of pack 2 wolves towards the RLUs}\nonumber\\
& &\mbox{made by pack 2.}\nonumber\\
\end{eqnarray*}
\begin{equation}
\tpa{u}=\nabla\cdot\left\{\ul{J}_{c_v}+\ul{J}_{d_v}+\ul{J}_{a_v}+\ul{J}_{b_v}
\right\}
\label{eq:w2a}
\end{equation}
where
\begin{displaymath}
\begin{array}{rclcc}
\ul{J}_{c_v}&=&v c_v(p) \frac{(\ul{x}-\ul{x}_v)}{|\ul{x}-\ul{x}_v|},&
c_v(0)\geq 0,&\zd{c_v}{p}\geq 0\\[2ex]
\ul{J}_{d_v}&=&d_v(v) \nabla v,&
d_v(0)\geq 0,&\zd{d_v}{v}\geq 0\\[2ex]
\ul{J}_{a_v}&=&va_v(p) \nabla p,&
a_v(0)\geq 0,&\zd{a_v}{p}\geq 0\\[2ex]
\ul{J}_{b_v}&=&-vb_v(q) \nabla q,&
b_v(0)\geq 0.& \\[2ex]
\end{array}
\end{displaymath}
Lastly, we require equations to describe the changes in RLU density ($p$ or $q$)
with time. The equation for $p$ is:
\begin{eqnarray*}
\lefteqn{\mbox{Rate of change in expected RLU density (pack 1)}}\nonumber\\
&=&
\mbox{Low level of continual marking}\nonumber\\
&+&
\mbox{Increased marking in the presence of RLUs from pack 2}\nonumber\\
&+&
\mbox{First order decay of RLUs with time}\nonumber
\end{eqnarray*}
\begin{equation}
\tpa{p}=u\left(l_p+m_pq\right)-f_pp,\label{eq:pa}\\[1ex]
\end{equation}
and the equation for $q$ is:
\begin{eqnarray}
\lefteqn{\mbox{Rate of change in expected RLU density (pack 2)}}\nonumber\\
&=&
\mbox{Low level of continual marking}\nonumber\\
&+&
\mbox{Increased marking in the presence of RLUs from pack 1}\nonumber\\
&+&
\mbox{First order decay of RLUs with time}\nonumber
\end{eqnarray}
\begin{equation}
\tpa{q}=v\left(l_q+m_qp\right)-f_qq.\label{eq:qa}
\end{equation}
Here $l_p$ and $l_q$ represent low level RLU
marking throughout the wolves' territory, $m_p$ and $m_q$ determine
the increase in RLU in the presence of competitive RLUs and $f_p$
and $f_q$ indicate decay dynamics for the RLUs, which are taken here to
be first order.
Mathematically we denote the study area by the
domain $\Omega$. A study of the system in two spatial dimensions
requires $\Omega\subset\RR^2$. Later in this paper we will mathematically
and numerically analyse the simpler one-dimensional case ($\Omega\subset\RR$).
To complete the mathematical formulation of our two-wolf pack model
we require boundary and initial
conditions. Biologically realistic boundary conditions may involve
local migration dynamics. The simplest possible boundary conditions
result, however, when we assume that wolves neither immigrate to nor emigrate
from the region $\Omega$. This situation is described by zero-flux boundary
conditions for $u$ and $v$, namely
\begin{equation}
\left\{\ul{J}_{c_u}+\ul{J}_{d_u}+\ul{J}_{a_u}+\ul{J}_{b_u}\right\}\cdot\ul{n}
\;\;\;\mbox{on $\partial\Omega$} \label{eq:w1_zf}
\end{equation}
and
\begin{equation}
\left\{\ul{J}_{c_v}+\ul{J}_{d_v}+\ul{J}_{a_v}+\ul{J}_{b_v}\right\}\cdot\ul{n}
\;\;\;\mbox{on $\partial\Omega$},\label{eq:w2_zf}
\end{equation}
where $\ul{n}$ is the outwardly oriented unit normal to the boundary of the
solution domain, $\partial\Omega$.
Initial conditions, describing the expected spatial distributions of wolves and
markings at the beginning of a study period, are given by
\begin{equation}
u(\ul{x},0)=u_0(\ul{x}),\;\;\;
v(\ul{x},0)=v_0(\ul{x}),\;\;\;
p(\ul{x},0)=p_0(\ul{x}),\;\;\;
q(\ul{x},0)=q_0(\ul{x}).
\label{eq:ics_1}
\end{equation}
We now show that the zero-flux boundary conditions
((\ref{eq:w1_zf})--(\ref{eq:w2_zf})) guarantee a constant
number of wolves for each pack within the domain $\Omega$.
At any given time, the total number of wolves from wolf pack 1 in
the domain $\Omega$ is
\begin{displaymath}
\int_\Omega u(\ul{x},t)\,d\ul{x}.
\end{displaymath}
The result follows directly from an application of the
divergence theorem. For wolf pack 1
\begin{eqnarray*}
\lefteqn{\frac{\partial}{\partial t}\int_\Omega u(\ul{x},t)\,d\ul{x}
=\int_\Omega \tpa{u}(\ul{x},t)\,d\ul{x}}\\
&=&\int_\Omega \nabla\cdot\left\{\ul{J}_{c_u}+\ul{J}_{d_u}+\ul{J}_{a_u}+\ul{J}_{b_u}\right\}\,d\ul{x}
=\int_{\partial\Omega} \left\{\ul{J}_{c_u}+\ul{J}_{d_u}+\ul{J}_{a_u}+\ul{J}_{b_u}\right\}\cdot\ul{n}\,ds
=0.
\end{eqnarray*}
An analogous argument holds true for pack 2.
Noting that the area of $\Omega$ is given by
\begin{displaymath}
A=\int_\Omega \,d\ul{x},
\end{displaymath}
we can calculate the average density of wolves from
pack 1 throughout the region $\Omega$ as
\begin{equation}
U_0=\frac 1 A\int_\Omega u_0(\ul{x})\,d\ul{x}
\label{eq:w1_av_dens}
\end{equation}
and the average density of wolves from pack 2 is
\begin{equation}
V_0=\frac 1 A\int_\Omega v_0(\ul{x})\,d\ul{x}.
\label{eq:w2_av_dens}
\end{equation}
To simplify further analysis we nondimensionalize the model system
(\ref{eq:w1a})--(\ref{eq:qa}) with boundary conditions
(\ref{eq:w1_zf})--(\ref{eq:w2_zf}) and intial conditions
(\ref{eq:ics_1})). This permits us to normalize the wolf density and domain
size as well as reduce the number of parameters. Defining $L=A^{1/n}$, where
$n$ is the dimension of the solution domain (either $n=1$ or $n=2$), we
write
\begin{displaymath}
u^*=\frac{u}{U_0},\;\;\;
v^*=\frac{v}{V_0},\;\;\;
p^*=\frac{pf_p}{U_0l_p},\;\;\;
q^*=\frac{qf_p}{V_0l_q},\;\;\;
t^*=tf_p,\;\;\;
\ul{x}^*=\frac{\ul{x}}{L}
\end{displaymath}
\begin{displaymath}
c_u^*=\frac{c_u}{Lf_p},\;\;\;
c_v^*=\frac{c_v}{Lf_p},\;\;\;
d_u^*=\frac{d_u}{L^2f_p},\;\;\;
d_v^*=\frac{d_v}{L^2f_p},\;\;\;
a_u^*=\frac{a_uU_0l_p}{L^2f_p^2},\;\;\;
a_v^*=\frac{a_vU_0l_p}{L^2f_p^2},
\end{displaymath}
\begin{displaymath}
b_u^*=\frac{b_uV_0l_q}{L^2f_p^2},\;\;\;
b_v^*=\frac{b_vV_0l_q}{L^2f_p^2},\;\;\;
\mu_p^*=\frac{m_pV_0l_q}{l_pf_p},\;\;\;
\mu_q^*=\frac{m_qU_0l_p}{l_qf_p},\;\;\;
\phi=\frac{f_q}{f_p}.
\end{displaymath}
For the nondimensionalized quantities to be well defined, we
implicitly assume that wolves from both packs
are present originally ($U_0>0$, $V_0>0$), that the domain $\Omega$ has a size
greater than zero ($L>0$), that both wolf packs have a non-zero low
level of RLU marking ($l_p>0$, $l_q>0$) and that the
RLU intensity decays with time ($f_p>0$).
Dropping the asterisks for notational simplicity, we write the
nondimensionalized system as (\ref{eq:w1a})--(\ref{eq:w2a}) and
\begin{equation}
\tpa{p}=u\left(1+\mu_p q\right)-p.\label{eq:p1a}
\end{equation}
\begin{equation}
\tpa{q}=v\left(1+\mu_q p\right)-\phi p.\label{eq:q1a}
\end{equation}
The boundary conditions (\ref{eq:w1_zf})--(\ref{eq:w2_zf}) remain unchanged,
and an appropriate nondimensionalization of the initial data,
\begin{displaymath}
u_0^*=\frac{u_0}{U_0},\;\;\;
v_0^*=\frac{v_0}{V_0},\;\;\;
p_0^*=\frac{p_0f_p}{U_0l_p},\;\;\;
q_0^*=\frac{q_0f_p}{V_0l_q},\;\;\;
\end{displaymath}
leaves the initial conditions (\ref{eq:ics_1}) unchanged
as well, once asterisks have been dropped.
Observe that our nondimensionalization
of space has rendered $\Omega$ an area ($\Omega\subset\RR^2$)
or length ($\Omega\subset\RR$) equal to unity.
Furthermore, with the nondimensionalization given here,
\begin{equation}
\int_{\Omega} u(\ul{x},t)\,d\ul{x}=
\int_{\Omega} v(\ul{x},t)\,d\ul{x}=1
\label{eq:int}
\end{equation}
and thus, at any given time, $u(\ul{x},t)$ and $v(\ul{x},t)$
are probability density functions for the location of wolves.
We believe that our quantitative,
behavior-based modelling exemplifies
an important new approach to understanding territorial
boundaries. In the next section we will show how
boundaries can result directly from solving
the two-wolf pack territorial model
They
depend intrinsically upon the behavioral parameters
and are thus not necessarily fixed.
This contrasts with previous models (see, for example,
Taylor and Pekins (1991))
which assume a fixed territorial boundary.
The difference becomes significant when questions, such as the stability
of the populations under stress, arise.
\subsection{Analytical and Numerical Results for the Territorial Model}
We now analyze the two-wolf pack territorial model
(\ref{eq:w1a})--(\ref{eq:w2a}),
(\ref{eq:w1_zf})--(\ref{eq:w2_zf}), (\ref{eq:ics_1}) and
(\ref{eq:p1a})--(\ref{eq:q1a}) mathematically and
numerically.
After making certain simplifying assumptions about the nonlinear
functions in (\ref{eq:w1a})--(\ref{eq:w2a}) we will calculate
steady state solutions in one spatial dimension. We then determine how
different sets of assumptions yield different qualitative forms for
the steady-state solutions.
One method of simplifying analysis is to consider cases where both
packs are similar and follow the same behavioral rules. This symmetry
between the two packs is then reflected in the nonlinear functions used in
(\ref{eq:w1a})--(\ref{eq:w2a}). In particular we will consider the case where
\begin{equation}
\begin{array}{ccc}
c_u(\cdot)=c_v(\cdot)=c(\cdot),&
d_u(\cdot)=d_v(\cdot)=d(\cdot),&
a_u(\cdot)=a_v(\cdot)=a(\cdot),\\
b_u(\cdot)=b_v(\cdot)=b(\cdot),&
\mu_p=\mu_q=\mu,&
\phi=1.
\end{array}
\label{eq:symm}
\end{equation}
The case of asymmetric pack interactions is very interesting from
both mathematical and biological perspectives and we anticipate
conducting further work on this.
As discussed earlier, our nondimensionalized space variable yields
a one-dimensional domain length equal to 1. The constraint (\ref{eq:int})
which was imposed by the nondimensionalization is thus re-written here in
one spatial dimension as
\begin{equation}
\int_0^1 u(x,t)\,dx=
\int_0^1 v(x,t)\,dx=1.
\label{eq:int_1D}
\end{equation}
As we are primarily interested
in the formation of territorial boundaries between packs, we assume that,
when dens (or organizing centers) are present, the dens for
each pack are at opposite ends of the domain. In particular, we choose
\begin{equation}
x_u=0,\;\;\;x_v=1.
\label{eq:org_cent}
\end{equation}
Under these assumptions we will show that:
\begin{enumerate}
\item When avoidance of competing RLUs is modelled
to include movement back towards the
den or organizing center ($\ul{J}_{c_u}$ and $\ul{J}_{c_v}$ are nonzero and
$\ul{J}_{a_u}$ and $\ul{J}_{a_v}$ are equal to zero
(see (\ref{eq:w1a})--(\ref{eq:w2a}))),
territorial pattern can arise with features that are qualitatively similar to
those observed in the field. In particular,
under the appropriate modelling assumptions,
low levels of expected wolf density and
high levels of expected RLU density are found close to the boundaries.
These correspond to the interpack `buffer zones' and to the
bowl-shaped RLU densities discussed elsewhere.
\item When the packs are modelled as having no den or organizing center,
and therefore avoidance of competing RLUs is modelled only by movement
down gradients of the intensity of the competing RLUs, patterns typically
{\em do not} arise; solutions are spatially constant.
(In this case, $\ul{J}_{c_u}$ and $\ul{J}_{c_v}$ are equal to zero and
$\ul{J}_{a_u}$ and $\ul{J}_{a_v}$ are nonzero.)
This result accentuates the importance of dens or
organizing centers for the generation of territorial patterns.
\end{enumerate}
\subsubsection{Time-Independent Solutions With Dens or Organizing Centers}
We now consider one-dimensional time-independent solutions to
(\ref{eq:w1a})--(\ref{eq:w2a}),
(\ref{eq:w1_zf})--(\ref{eq:w2_zf}), (\ref{eq:ics_1}) and
(\ref{eq:p1a})--(\ref{eq:q1a}), (\ref{eq:symm})--(\ref{eq:org_cent}),
with
\begin{equation}
a(u)=a(v)=0,
\label{eq:a_zero}
\end{equation}
so that $\ul{J}_{a_u}$ and $\ul{J}_{a_v}$ are equal to zero.
Thus the avoidance of competitive RLUs is solely through the directed
movement (or `convection') back to the den at a speed given by
$c(q)$ (pack 1) or $c(p)$ (pack 2).
For the purpose of mathematically analyzing the solutions we
choose the simplest possible form for $c(\cdot)$; that is, we choose
\begin{equation}
c(q)=cq,\;\;\; c(p)=cp\;\;\;\mbox{(c constant)}.
\label{eq:cdef}
\end{equation}
Initially we consider the case where RLU marking is independent of the
level of competing RLUs ($\mu=0$). We then relax this restriction
and consider how patterns can evolve for nonzero $\mu$.
\begin{center}
\large{\underline{The case $\mu=0$}}
\end{center}
When $\mu=0$ the time-independent solutions to (\ref{eq:p1a})--(\ref{eq:q1a})
are simply given by $p=u$ and $v=q$. Thus steady state solutions to
(\ref{eq:w1a})--(\ref{eq:w2a}), (\ref{eq:symm})--(\ref{eq:cdef}) satisfy
\begin{eqnarray*}
0&=&\pax\left\{cuv+d(u)u_x-ub(u)u_x\right\}\\
0&=&\pax\left\{-cuv+d(v)v_x-vb(v)v_x\right\}.
\end{eqnarray*}
Application of the zero-flux boundary conditions
(\ref{eq:w1_zf})--(\ref{eq:w2_zf}) yields
\begin{eqnarray}
0&=&cuv+d(u)u_x-ub(u)u_x\label{eq:u_mu_0}\\
0&=&-cuv+d(v)v_x-vb(v)v_x.\label{eq:v_mu_0}
\end{eqnarray}
Note that these equations are symmetric under the transformation
$x\rightarrow 1-x$, $u\leftrightarrow v$. This symmetry, and the constraint
(\ref{eq:int_1D}), leads us to expect that the time-independent solutions,
$u(x)$ and $v(x)$, are mirror images of each other about the mid-point
$x=1/2$, so that $u(0)=v(1)$ and $u(1)=v(0)$.
We now simplify (\ref{eq:u_mu_0})--(\ref{eq:v_mu_0}) for the case where
the diffusive flux dominates the aggregative flux that arises
through following familiar RLUs; that is, $d(w)>wb(w)$ for all $w\geq 0$.
Defining
\begin{equation}
\Gamma(w)=\int_0^w\left\{d(\omega)-\omega b(\omega)\right\}\,d\omega,
\label{eq:gam_def}
\end{equation}
we note that $\Gamma(w)$ is a non-negative invertible function for
$w\geq 0$ because $\Gamma(0)=0$ and $\Gamma '(w)=d(w)-wb(w)>0$.
Adding (\ref{eq:u_mu_0}) and (\ref{eq:v_mu_0}) and integrating
with respect to $x$ we observe that
\begin{equation}
\Gamma(u)+\Gamma(v)=\Gamma(u(0))+\Gamma(v(0))=k\;\;\;\mbox{(constant)}
\label{eq:gam_cons}
\end{equation}
and thus (\ref{eq:gam_cons}) is conserved for all $u(x)$ and $v(x)$
satisfying (\ref{eq:u_mu_0})--(\ref{eq:v_mu_0}). We use the fact that
$\Gamma$ is invertible to separate $u$ and $v$ in
(\ref{eq:u_mu_0})--(\ref{eq:v_mu_0}) and rewrite them as
\begin{eqnarray}
0&=&cu\Gamma^{-1}(k-\Gamma(u))+(d(u)-ub(u))u_x\label{eq:u_mu_0a}\\
0&=&-cv\Gamma^{-1}(k-\Gamma(v))+(d(v)-vb(v))v_x\label{eq:v_mu_0a}.
\end{eqnarray}
Integration of (\ref{eq:u_mu_0a}) yields $u(x)$ and similarly
integration of (\ref{eq:v_mu_0a}) yields $v(x)$.
By way of examples, we now consider two special cases:
\begin{enumerate}
\item \underline{$d(w)=d$, (constant) and $b(w)=0$:}
In is a case, the wolves move randomly (via simple diffusion) in the
absence of a competing pack.
Here equation (\ref{eq:gam_def})
yields $\Gamma(w)=dw$ and thus $\Gamma^{-1}(w)=d^{-1}w$.
Therefore equation (\ref{eq:gam_cons})
indicates that $u+v=k/d$ is conserved across the solution domain $0\leq x \leq 1$.
Equations (\ref{eq:u_mu_0a}) and (\ref{eq:v_mu_0a}) become
\begin{eqnarray}
u_x&=&\frac{-c}{d}u\left(\frac{k}{d}-u\right)\label{eq:log_u}\\
v_x&=&\frac{c}{d}v\left(\frac{k}{d}-v\right).\label{eq:log_v}
\end{eqnarray}
These are simply versions of the logistic equations with
space as the independent variable. Explicit solutions for
$u$ and $v$ are given by
\begin{eqnarray}
u(x)=\frac{ku(0)\exp\{-\frac{d^2}{ck}x\}}{k+du(0)[\exp\{-\frac{d^2}{ck}x\}-1]}
\label{eq:log_u_sol}\\
v(x)=\frac{kv(0)\exp\{\frac{d^2}{ck}x\}}{k+dv(0)[\exp\{\frac{d^2}{ck}x\}-1]}
\label{eq:log_v_sol}
\end{eqnarray}
where $k$ is $k=d(u(0)+v(0))$ (see (\ref{eq:gam_cons})).
By integrating (\ref{eq:log_u_sol}) and (\ref{eq:log_v_sol}) with respect
to $x$ and applying the constraint (\ref{eq:int_1D})
we have
\begin{eqnarray}
u(0)&=&(u(0)+v(0)) \frac{\exp\left\{-\frac{d}{c(u(0)+v(0))^2}\right\}-1}
{\exp\left\{-\frac{d}{c(u(0)+v(0))}\right\}-1}\label{eq:u0_int1}\\
v(0)&=&(u(0)+v(0)) \frac{\exp\left\{\frac{d}{c(u(0)+v(0))^2}\right\}-1}
{\exp\left\{\frac{d}{c(u(0)+v(0))}\right\}-1}.\label{eq:v0_int1}
\end{eqnarray}
These two equations can be solved simultaneously to yield the two
unknown quantities, $u(0)$ and $v(0)$, in terms of the parameters $c$ and $d$.
Solutions are then given by (\ref{eq:log_u_sol})--(\ref{eq:log_v_sol})
where $k=d(u(0)+v(0))$.
A numerical solution to the related time-dependent problem
((\ref{eq:w1a})--(\ref{eq:w2a}),
(\ref{eq:w1_zf})--(\ref{eq:w2_zf}), (\ref{eq:ics_1}) and
(\ref{eq:p1a})--(\ref{eq:q1a}), (\ref{eq:symm}),
(\ref{eq:org_cent}), (\ref{eq:cdef}), $a(w)=0$, $d(w)=d$, $b(w)=0$,
and $\mu=0$) achieved the predicted heterogeneous steady-state pattern
for large time (Figure 2).
Initial conditions were given by
\begin{equation}
u(x,0)=1-H(1/2),\;\;\;v(x,0)=H(1/2).
\label{eq:H_ic}
\end{equation}
This numerical solution confirms our analytical
prediction that $u+v$ is constant for all $x$. We thus observe that
distinct territorial patterns may form. However this solution differs
qualitatively from observed results in that
\begin{enumerate}
\item The total expected wolf density
($u+v$) remains constant throughout and is not reduced near the inter-pack
boundary and
\item The expected density of RLU markings ($p=u$ and $q=v$) is
no higher near territorial boundaries than near the den. We will show later
in this section that this deficiency in the results is rectified when
$\mu\neq 0$.
\end{enumerate}
Before considering the $\mu\neq 0$ case, we briefly analyse
a case with $\mu=0$ and nonlinear diffusion.
\item \underline{$d(w)=dw$ ($d$ constant), and $b(w)=b$ ($b$ constant), $d>b$:}
Here equation (\ref{eq:gam_def}) yields
$\Gamma(w)=(d-b)w^2/2$ and thus $\Gamma^{-1}(w)=\sqrt{2w/(d-b)}$.
Therefore equation (\ref{eq:gam_cons})
says that $u^2+v^2=2k/(d-b)$ is conserved across the solution domain
$0\leq x \leq 1$.
Equations (\ref{eq:u_mu_0a}) and (\ref{eq:v_mu_0a}) become
\begin{eqnarray}
0&=&u\left[c\left(\frac{2k}{d-b}-u^2\right)^{1/2}+(d-b)u_x\right]\label{eq:asin_u}\\
0&=&v\left[-c\left(\frac{2k}{d-b}-v^2\right)^{1/2}+(d-b)v_x\right].\label{eq:asin_v}
\end{eqnarray}
Integration of (\ref{eq:asin_u})--(\ref{eq:asin_v}) and application of the
integral constraint (\ref{eq:int_1D}) yield the following solution
which is valid when $d-b\geq 2c/\pi$:
\begin{eqnarray}
u(x)=\left\{
\begin{array}{ll}
\frac{2}{1+(d-b)(2-\pi/2)/c}&0\leq x\leq x_0\\[3ex]
\sqrt{\frac{2k}{d-b}}\cos\left(\frac{c(x-x_0)}{d-b}\right)&x_0 0$. Thus steady state solutions to
(\ref{eq:w1a})--(\ref{eq:w2a}), (\ref{eq:symm})--(\ref{eq:cdef}) satisfy
\begin{eqnarray*}
0&=&\pax\left\{cuq+d(u)u_x-ub(p)p_x\right\}\\
0&=&\pax\left\{-cvp+d(v)v_x-vb(q)q_x\right\}
\end{eqnarray*}
with $p$, $q$, $p_x$ and $q_x$ defined by (\ref{eq:p_ss})--(\ref{eq:qd_ss}).
Application of the zero-flux boundary conditions
(\ref{eq:w1_zf})--(\ref{eq:w2_zf}) yields
\begin{eqnarray}
0&=&cuq+d(u)u_x-ub(p)p_x\label{eq:u_mu_n0}\\
0&=&-cup+d(v)v_x-vb(q)q_x\label{eq:v_mu_n0}
\end{eqnarray}
Thus (\ref{eq:p_ss})--(\ref{eq:v_mu_n0}) defines a system of two
coupled nonlinear ordinary differential equations for $u$ and $v$. Therefore,
given any prescribed $d(\cdot)$ and $b(\cdot)$ it would be possible
to analyse the behavior of this system in the $u$-$v$ phase plane.
Our approach is to mathematically analyse the simplest possible case in depth,
$d(w)=d$ (constant) and $b(w)=0$. We anticipate that
the behaviors exhibited in this simple case may also
be present for more complex $d(\cdot)$ and $b(\cdot)$.
Substitution of (\ref{eq:p_ss}) and (\ref{eq:q_ss}) into
(\ref{eq:u_mu_n0})--(\ref{eq:v_mu_n0}) yields
\begin{eqnarray}
0&=&c\frac{uv(1+\mu u)}{1-\mu^2uv}+du_x\label{eq:u_fun}\\
0&=&-c\frac{uv(1+\mu v)}{1-\mu^2uv}+dv_x\label{eq:v_fun}
\end{eqnarray}
We can now formulate a conservation law
as we did earlier for the $\mu=0$ case (see (\ref{eq:gam_cons})).
Defining
\begin{displaymath}
\Gamma_1(w)=\int_0^w \frac {d\omega}{1+\mu\omega}=\log(1+\mu w)/\mu
\end{displaymath}
we observe from (\ref{eq:u_fun})--(\ref{eq:v_fun}) that
\begin{displaymath}
\Gamma_1(u)+\Gamma_1(v)=\Gamma_1(u(0))+\Gamma_1(v(0))=k\;\;\;\mbox{(constant)}.
\end{displaymath}
Thus
\begin{equation}
(1+\mu u)(1+\mu v)=\exp(\mu k)
\label{eq:gam1_cons}
\end{equation}
describes $u$ in terms of $v$ and vica-versa.
Substitution of (\ref{eq:gam1_cons}) into (\ref{eq:v_fun}) yields
\begin{equation}
0=-c\frac{v\left(\frac{\exp(\mu k)}{1+\mu v} -1\right)/\mu}
{1-\mu v\left(\frac{\exp(\mu k)}{1+\mu v} -1\right)}
+\frac{dv_x}{1+\mu v}
\label{eq:v_int}
\end{equation}
while substitution into (\ref{eq:u_fun}) yields a similar equation for $u$.
Equation (\ref{eq:v_int}) is separable, and integration yields
an implicit formula for $v(x)$, namely
\begin{displaymath}
\frac{v(x)}{(E-\mu v(x))(1+\mu v(x))^E}=\frac{v(0)\exp\{\frac{cEx}{d\mu}\}}
{(E-\mu v(0))(1+\mu v(0))^E},
\end{displaymath}
where we have written $E=\exp(\mu k) -1$.
For the special case $E=1$ we can calculate the solution explicitly as
\begin{eqnarray}
v(x)&=&\frac 1 {2\mu^2}\left[\left(\mu^2v(0\right)^2-1)\exp\left\{\frac{-cx}{d\mu}\right\}/v(0)\right.
\nonumber\\
&+&\left.\left(\left(\mu^2v(0)^2-1\right)^2\exp\left\{\frac{-2cx}{d\mu}\right\}/v(0)^2+4\mu^2\right)^{1/2}
\right].\label{eq:v_expl}
\end{eqnarray}
Noting the symmetry of (\ref{eq:u_fun})--(\ref{eq:v_fun}) under the
transformation $x\rightarrow 1-x$ and $u\leftrightarrow v$ we can deduce
that $u(0)=v(1)$. From equation (\ref{eq:gam1_cons}) we observe that,
when $E=1$, $(1+\mu u(0))(1+\mu v(0))=2$, and thus
\begin{equation}
(1+\mu v(1))(1+\mu v(0))=2.\label{eq:Newt1_it}
\end{equation}
Integration of (\ref{eq:v_expl}) with respect to $x$
and application of the constraint (\ref{eq:int_1D}) defines $v(0)$ in
terms of $c/d$, and $\mu$. Specifically, $v(0)$ must satisfy
\begin{eqnarray}
1&=&-\frac{dv(0)}{2\mu c}\left[1-\exp\left\{\frac{-c}{d\mu}\right\}\right]
+\frac d c \left[-
\left(\left(\frac{v(0)}{2\mu}\right)^2\exp\left\{\frac{-2c}{d\mu} \right\}+1\right)^{1/2}\right.\nonumber\\
& &-\frac 1 2\log\left(1-
\left(\left(\frac{v(0)}{2\mu}\right)^2\exp\left\{\frac{-2c}{d\mu} \right\}+1\right)^{-1/2}
\right)\nonumber\\
& & +\frac 1 2\log\left(1+
\left(\left(\frac{v(0)}{2\mu}\right)^2\exp\left\{\frac{-2c}{d\mu} \right\}+1\right)^{-1/2}
\right)\nonumber\\
%\end{eqnarray*}
%\begin{eqnarray*}
& &+\left(\left(\frac{v(0)}{2\mu}\right)^2+1\right)^{1/2}
+\frac 1 2\log\left(1-
\left(\left(\frac{v(0)}{2\mu}\right)^2+1\right)^{-1/2}
\right)\nonumber\\
& &\left.-\frac 1 2\log\left(1+
\left(\left(\frac{v(0)}{2\mu}\right)^2+1\right)^{-1/2}
\right) \right].\label{eq:Newt2_it}
\end{eqnarray}
We have thus derived a one-parameter family of solutions to
(\ref{eq:u_fun})--(\ref{eq:v_fun}) for the case $E=1$;
when given a value for $\mu$, we can calculate $v(0)$ and $c/d$ by simultaneously
solving two equations, the first of which is
equation (\ref{eq:Newt1_it}) in conjunction with (\ref{eq:v_expl})
for the definition of $v(1)$ in terms of $v(0)$, and the second of
which is equation (\ref{eq:Newt2_it}).
In Table 1 we computed the values of $v(0)$ and
$c/d$, while varying $\mu$. The corresponding profiles for $v(x)$,
as given by (\ref{eq:v_expl}), are shown in Figure 4. Finally, the
cumulative density for both wolf-packs and for both sets of RLUs are
shown in Figure 4. Note that the cumulative wolf-pack density
drops near the boundary at $x=1/2$ and the cumulative RLU density
is highest there. Thus this solution reflects the biological observations
(discussed elsewhere) more accurately than the $\mu=0$ case.
A numerical solution of the related time-dependent problem
((\ref{eq:w1a})--(\ref{eq:w2a}),
(\ref{eq:w1_zf})--(\ref{eq:w2_zf}), (\ref{eq:ics_1}) and
(\ref{eq:p1a})--(\ref{eq:q1a}), (\ref{eq:symm}),
(\ref{eq:org_cent}), (\ref{eq:cdef}), $a(w)=0$ and $b(w)=0$)
for the $\mu=0.48$ case from Table 1
($c=1$ and $d=1/7.3584$) achieved the steady state pattern for
large time (Figure 6). Initial conditions were given by (\ref{eq:H_ic}).
(Compare the numerically calculated values for
$u$ (Figure 6) with the analytically predicted pattern (Figure 4, top
curve).)
Note that, even though the cumulative RLU
densities have a maximum at $x=1/2$ (Figure 5),
the individual RLU densities {\em do not} have such a maximum (Figure 6).
We can show that this must indeed be the
case by looking at equations (\ref{eq:pd_ss}) and (\ref{eq:qd_ss}) in
conjuction with equation (\ref{eq:gam1_cons}). By differentiating
(\ref{eq:gam1_cons}) with respect to $x$ we observe that
\begin{displaymath}
0=\frac{v_x}{1+\mu v}+\frac{u_x}{1+\mu u}
\end{displaymath}
and thus (\ref{eq:pd_ss}) and (\ref{eq:qd_ss}) can be simplified to give
\begin{eqnarray}
p_x&=&\frac{(1+\mu u)(\mu u -1)}{(1-\mu u v)^2} v_x\label{eq:pd_sim}\\
q_x&=&\frac{(1+\mu v)(\mu v -1)}{(1-\mu u v)^2} u_x.\label{eq:qd_sim}
\end{eqnarray}
Thus we can calculate the critical points for $p$ and $q$ as
\begin{eqnarray}
p_x=0 \Leftrightarrow v_x=0\;\;\;\mbox{or}\;\;\;u=1/\mu\label{eq:pxz}\\
q_x=0 \Leftrightarrow u_x=0\;\;\;\mbox{or}\;\;\;v=1/\mu.\label{eq:qxz}
\end{eqnarray}
Equations (\ref{eq:u_fun}) and (\ref{eq:v_fun}) indicate that the
signs of $u_x$ and $v_x$ do not change. (You may recall our earlier assumption
that $1-\mu^2 uv>0$). Thus an interior maximum for $p$ requires that
$u=1/\mu$ and an interior maximum for $q$ requires that $v=1/\mu$.
Substitution of $u=1/\mu$ and $E=exp(\mu k) -1$
into (\ref{eq:gam1_cons}) indicates that an interior maximum
for $p$ can only occur if
\begin{equation}
u=\frac 1 u,\;\;\;v=\frac{E-1}{2\mu}
\label{eq:pmax}
\end{equation}
and substitution of $v=1/\mu$ and $E=exp(\mu k) -1$
into (\ref{eq:gam1_cons}) indicates that an interior maximum
for $q$ can only occur if
\begin{equation}
u=\frac{E-1}{2\mu},\;\;\;v=\frac 1 u.
\label{eq:qmax}
\end{equation}
Thus for the case $E=1$, an interior maximum for $p$ requires that
$v=0$, however never occurs (see equation (\ref{eq:v_expl})).
We thus conclude that there are no interior
maximums for $p$ when $E=1$. A similar argument applies for $q$.
Based on the above argument and equations (\ref{eq:pmax}) and
(\ref{eq:qmax}), it is natural to ask whether interior maximums
for $p$ and $q$ are possible when $E>1$.
If $p$ were to assume an interior maximum, then it's value would be given
by (\ref{eq:p_ss}) evaluated at (\ref{eq:pmax}), namely
\begin{equation}
p_{\mbox{max}}=\frac{1+E}{\mu(3-E)}.
\end{equation}
Similarly an interior maximum for $q$ would be given by (\ref{eq:q_ss})
evaluated at (\ref{eq:qmax}), namely $q_{\mbox{max}}=p_{\mbox{max}}$.
If bounded interior maximums $p$ and $q$ were to occur, they would thus be for
$E$ in the range $1