\xpa{u}=-\xpa{v},\;\;\xpa{p}=-\xpa{q}, \end{displaymath} \begin{displaymath} \pax(u+v)=0,\;\;\xpa{q}=\frac{\xpa{v}(1+m(p))}{1+vm'(p)}>0. \end{displaymath} and \begin{eqnarray*} \lefteqn{(u+v)_{xx}}\\ &=&\frac 1 d \left\{c(p)v-c(q)u\right\}_x\\ &=& \frac 2 d\left\{c'(p)up_x-c(p)u_x\right\}\\ &=& \frac {2u^2} {d}\pax\left\{\frac{c(p)}{p}\frac{p}{u}\right\}= \frac {2u^2} {d}\pax\left\{\frac{c(p)}{p}(1+m(q))\right\}\\ &=& \frac {2u^2} {d}\left\{\paz{p}\left(\frac{c(p)}{p}\right) (1+m(q))p_x+\frac{c(p)}{p}m'(q)q_x\right\}. \end{eqnarray*} A sufficient condition for the right hand side of the last line to be positive is that $c(p)$ is convex. In this case $x=1/2$ is a minimum for $u+v$ corresponding to a buffer zone for the interacting packs. \section{Dependence of Territories on Behavioral Responses} In this section we consider the qualitative dependence of the wolf territories on the form of the movement response function $c$ and the scent-marking response function $m$. Functional responses have been widely analysed in the context of predator-prey dynamics \cite{Holling:MESC-66-1}. The forms we use here for the behavioral responses are analogous to the Holling type II and III functional responses; the former is a convex monotonic function while the latter exhibits switching. We continue to focus on the model describing two identical interacting pack, with dens at opposing ends of the domain and assume no explicit spatial dependence in the movement response function $c$ (\ref{eq:su})--(\ref{eq:sq}). Using piecewise linear caricatures for the functions we ask how the wolf territories and scent-marking densities depend on the behavioral response. \subsection{Case I: No Marking Response to Foreign RLUs} This case is characterized by no scent marking response of wolves to foreign RLUs ($m=0$) and a linear increase in the directed movement rate in response to foreign RLUs ($c(q)=\gamma q$) (Figure 4(a)). It is assumed that the function $c$ is none-the-less bounded, but that $q,p\leq c_\infty/\gamma$ (Figure 4(a)). Here (\ref{eq:sp}) and (\ref{eq:sq}) simply yield $p=u$ and $q=v$. Substitution into (\ref{eq:su}) and (\ref{eq:sv}), and addition of these equations subject to the constraint (\ref{eq:icon}) yields $u(x)+v(x)=2$ pointwise on $0\leq x \leq 1$. Thus the territories are given by solutions to two logistic equations with space as the independent variable: \begin{equation} u_x=\frac{-\gamma}{d}u\left(2-u\right),\;\; v_x=\frac{\gamma}{d}v\left(2-v\right) \label{eq:logistic} \end{equation} (Figure 5). The initial conditions $u(0)$ and $v(0)$ to this system are chosen to satisfy the integral constraints (\ref{eq:icon}). We can construct an energy function on the globally stable invariant manifold $p=u$, $q=v$: \begin{eqnarray} F&=&\int_0^1\phi^2\,dx,\;\;\phi=u+v-2,\label{eq:F}\\ \tpa{F}&=&-2\int_0^1{\xpa{\phi}}^2\,dx\label{eq:dFdt}. \end{eqnarray} Thus the territory corresponds to a state of lowest potential, where the potential is measured as the squared deviation of away from a constant cumulative wolf density $u+v=2$. The function $F$ is not strictly a Lyapunov function because it is semi-definite and is restricted to the stable subspace $p=u$, $q=v$. The territories calculated in this section (Figure 5) differ qualitatively from results observed in the field 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. \end{enumerate} We now consider a modification of the above case which rectifies these two discrepancies. \subsection{Case II: Marking Response to Foreign RLUs} Here we assume that $m(q)$ and $c(q)$ have similar functional forms (Figures 4(a) and 4(b)). The part of the functions of interest are given by $m(q)=\mu q$ and $c(q)=\gamma q$. Although the functions $m$ and $c$ in Figures 4(a) and 4(b) are bounded, we assume that that $q,p\leq c_\infty/\gamma$ and $q,p\leq m_\infty/\mu$ (Figures 4(a) and 4(b)). Here (\ref{eq:sp}) and (\ref{eq:sq}) yield \begin{eqnarray} p&=&\frac{u(1+\mu v)}{1-\mu^2uv}\label{eq:p_ss}\\ q&=&\frac{v(1+\mu u)}{1-\mu^2uv}\label{eq:q_ss} \end{eqnarray} and substitution into (\ref{eq:su}) and (\ref{eq:sv}) yields \begin{eqnarray} 0&=&\gamma\frac{uv(1+\mu u)}{1-\mu^2uv}+du_x\label{eq:u_fun}\\ 0&=&-\gamma\frac{uv(1+\mu v)}{1-\mu^2uv}+dv_x\label{eq:v_fun}. \end{eqnarray} Defining \begin{displaymath} \Gamma(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{equation} \Gamma(u)+\Gamma(v)=\Gamma(u(0))+\Gamma(v(0))=k(u(0),v(0))\;\;\;\mbox{(constant)}. \label{eq:constant} \end{equation} 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 vice-versa. Substitution of (\ref{eq:gam1_cons}) into (\ref{eq:v_fun}) yields \begin{equation} 0=-\gamma\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{equation} \frac{v(x)}{(E-\mu v(x))(1+\mu v(x))^E}=\frac{v(0)\exp\{\frac{\gamma Ex}{d\mu}\}} {(E-\mu v(0))(1+\mu v(0))^E}, \label{eq:E} \end{equation} where we have written $E=\exp(\mu k(u(0),v(0)) -1$, and the similar $u$ equation can be solved by the same method. For the special case $E=1$ the solution can be calculated explicitly. Finally, the initial conditions $u(0)$ and $v(0)$ must be chosen so as to satisfy the integral constraints (\ref{eq:icon}). Profiles for the pack and scent mark densities are shown in Figure 6(a) and the cumulative expected densities for both wolf packs and for both scent marks are shown in Figure 6(b). Note that the cumulative wolf-pack density drops near the territorial boundary at $x=1/2$ and the cumulative RLU density is highest there. Thus this solution reflects the biological observations more accurately than the $m=0$ case. We use equation (\ref{eq:gam1_cons}) to observe that \begin{displaymath} 0=\frac{v_x}{1+\mu v}+\frac{u_x}{1+\mu u}. \end{displaymath} This can be used to simplify expressions for $p_x$ derived from (\ref{eq:p_ss}) and $q_x$ derived from (\ref{eq:q_ss}): \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*} Because $u(x)$ and $v(x)$ have no critical points (Section 4) possible interior maximums for $p(x)$ and $q(x)$ are given when $u(x)=1/\mu$ and $v(x)=1/\mu$, respectively. Thus there is an interior maximum for $p$ if and only if $u(0)\geq 1/\mu\geq u(1)$ and there is an interior maximum for $q(x)$ if and only if $v(0)\leq 1/\mu\leq v(1)$ (Figure 6(a)). In other words, if the behavioral response function $m$ is sufficiently steep then $1/\mu$ will lie in the above interval and bowl-shaped scent marking densities will arise. \subsection{Case III: Switching in Movement Response to Foreign RLUs} We now include `switching' in the movement response of wolves to foreign scent marks. Specifically, we assume that there is no response at all until the foreign scent mark has reached a critical value $q_c$ at the point $x_c$, after which there is a return to the den at speed $c_\infty$. The marking response is chosen to be the same as in Case II above. Hence we have $m(q)=\mu q$ (Figure 4(b)) and $c(q)=c_\infty H(q-q_c)$ where $H(\cdot)$ is the Heaviside function (Figure 4(c)). Using symmetry, equations (\ref{eq:su}) and (\ref{eq:sv}) are solved to give \begin{eqnarray} u&=&\left\{ \begin{array}{ccc} u(0)&\mbox{if}&0\leq x \leq x_c\\ u(0)\exp\left(-\frac{c_\infty}{d} (x-x_c)\right)&\mbox{if}&x_c< x\leq 1 \end{array} \right.\label{eq:usw}\\ v&=&\left\{ \begin{array}{ccc} u(0)&\mbox{if}&1-x_c\leq x \leq 1\\ u(0)\exp\left(-\frac{c_\infty}{d} (1-x-x_c)\right)&\mbox{if}&0\leq x < 1-x_c \end{array} \right.\label{eq:vsw} \end{eqnarray} where $u(0)$ is given in terms of $x_c$ by the integral constraint (\ref{eq:icon}) as \begin{equation} u(0)\left(x_c+\frac{d}{c_\infty}\left(1-\exp\left(\frac{c_\infty}{d} (x_c-1)\right)\right)\right)=1. \label{eq:icon1} \end{equation} The equations for $p$ and $q$ are given by (\ref{eq:p_ss}) and (\ref{eq:q_ss}). Using (\ref{eq:q_ss}) and (\ref{eq:usw})--(\ref{eq:vsw}) we calculate \begin{eqnarray} q_c&=&q(x_c)=\frac{v(x_c)(1+\mu u(0))}{1-\mu^2u(0)v(x_c)}\label{eq:q_c}\\ v(x_c)&=& u(0)\exp\left(-\frac{c_\infty}{d}(1-2x_c)\right).\nonumber \end{eqnarray} Simultaneous solution of (\ref{eq:icon1}) and (\ref{eq:q_c}) yields the maximum expected density of wolves from a pack ($u(0)$) and the width of the buffer zone ($1-2x_c$) in terms of the parameters $c_\infty$, $d$, $q_c$ and $\mu$. Whereas some buffer zone is evident without switching (Figure 6(b)), the switching function can greatly enhance the the magnitude of the buffer zone (Figure 7). For any given width, the depth of the buffer zone is an increasing function of the ratio of directed to random motion in the presence of foreign scent marks $c_\infty/d$. \subsection{Case IV: Switching in Marking Response to Foreign RLUs} Lastly we consider the case where there is switching in the scent marking behavioral response. The movement response is chosen to be the same as in Case II above. The functional forms are thus $m(q)=m_\infty H(q-q_m)$ (Figure 4(d)) and $c(q)=\gamma q$ (Figure 4(a)). Figure 8 shows the stable time-independent equilibria for (\ref{eq:p_ss})--(\ref{eq:q_ss}) in terms of $u$ and $v$. Note that if $(u,v)$ lies within the shaded region of Figure 8 there is not a unique time-independent solution. However, providing this region in $u-v$ phase space is bypassed the solution can be calculated in each sub-region in a manner similar to that used in Case II above. At the transition from Region B to Region D the expected scent mark density for pack 1 ($p(x_m)$) jumps up from $u(x_m)$ to $(1+m_\infty)u(x_m)$ and at the transition from Region D to Region C the expected scent mark density for pack 2 ($q(1-x_m)$) jumps down from $(1+m_\infty)v(1-x_m)$ to $v(1-x_m)$ (Figure 8). Using symmetry, the relationship between $u(x)$ and $v(x)$ given in terms of $k=u(0)+(1+m_\infty)v(0)-m_\infty q_m$ is \begin{equation} \begin{array}{ccccc} u+v+m_\infty v&=& k+m_\infty q_m& \mbox{if}&0\leq x \leq x_m\\ u+v&=& k&\mbox{if}&x_m< x < 1-x_m\\ u+v+m_\infty u&=& k+m_\infty q_m &\mbox{if}&1-x_m\leq x \leq 1\\ \end{array} \label{eq:mswitch} \end{equation} The value of $k$ is chosen so as to satisfy the integral constraints (\ref{eq:icon}). A sample trajectory is shown in Figure 8 and the corresponding solution is shown in Figure 9. Here the step function form for $m(q)$ yields a solution which is discontinuous in $p$ and $q$. In other words, an abrupt jump in $p$ and $q$ at the edge of the buffer zone corresponding to a sharp lip on the bowl of scent marking densities is direct consequence of switching in $m$. In summary, our analysis of the piecewise linear system response functions has demonstrated how the functions determine the qualitative nature of the territories. Specifically \begin{enumerate} \item Territories can form with a convex movement response function and no scent marking response function $(m=0)$ but they typically do not exhibit the bowl shaped scent mark densities and the buffer zones between territories. \item The addition of a convex scent marking response function can give rise to the bowl shaped scent marking densities and buffer zones. A condition for the bowl shaped scent marking densities is that the scent marking response function $(m)$ is sufficiently steep. \item Inclusion of switching in the movement response function can give rise to very distinct, wide buffer zones. \item Inclusion of switching in the scent marking response function can give rise to very distinct, sharp edges on the scent mark bowl. \end{enumerate} \section{Numerical Solution in Two Spatial Dimensions} The full nonlinear system (\ref{eq:w1a})--(\ref{eq:qa}) was solved numerically using finite differences in two spatial dimensions. Computation was terminated when the solution had effectively reached a stationary profile. This stationary profile corresponds to territories, as described above. Shown are the two- and three-pack interactions (Figure 10). The numerical solution over a range of smooth movement behavior functions ($c$) and scent marking behavior functions ($m$) indicated that switching in $c$ is needed for distinct territories and switching in $m$ is needed for an lip on the bowl of scent marking densities. When there was no switching in $m$ the regions of high scent marking were not localized to the edges of the `bowl' but formed large two-dimensional hot-spots in a `blotchy' manner. The precise forms of the switching functions was not crucial as long as they were given qualitatively by Holling type III. This contrasts with the 1 dimensional results illustrated in Figure 6; there the bowl-shaped RLU densities arise with a convex $m$ and switching in $m$ only serves to sharpen the lip. \section{Discussion} In this paper we have analysed a simple mechanistic mathematical model which assumes that wolf movement and scent marking is mediated by the presence or absence of foreign scent marks. Specifically, upon encountering a foreign scent mark a wolf is assumed to increase its own scent marking and move toward the organizational center of its pack. The steady states that result correspond to territories. We have shown that, under a wide variety of assumptions about the movement and scent marking response functions, territories are determined uniquely and buffer zones exist between territories. The presence or absence of bowl shaped scent mark densities observed in nature specifically depends upon the shape and steepness of the scent marking response function. Numerical experimentation in two spatial dimensions indicates that switching in both the movement response function and in the scent marking response function are needed to generate realistic territories with buffer zones and bowl shaped scent mark densities. Although piecewise constant switching functions were used for the mathematical analysis, the numerical solutions with smooth switching functions indicate that it is the qualitative form of the function (switching vs.\ no switching) that determines the kind of territorial patterns. While analysis was primarily restricted to the explicit calculation of territories as time-independent solutions, all of our numerical results indicate that the territories are stable configurations. We conjecture that the territories are globally stable when unique, and have local stability properties when switching in $m$ causes loss of uniqueness (see Figures 2(b) and 8). Analytical approaches to stability would be a worthwhile topic for additional research. These could possibly use an energy method similar to that in Section 6.1. Although nonlinear partial differential equations have been applied successfully to ecological problems \cite{Holmes:EC-94-17,Okubo:DEP80} including the spatial segregation of interacting species \cite{Shigesada:JTB-79-83,Mimura:1980:SSC} and understanding social aggregations \cite{Gruenbaum:1994:MSA} the modelling given here represents a new approach to understanding behavioral aspects of territory formation. A significant feature of our study is that the seemingly complex formation of wolf territories can be reduced to a relatively simple set of formulas involving scent marking.