1 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 expected densities of wolves at point x and time t, as direct field observations typically will not yield the exact densities. Thus, the pertinent state variables are:

u(x, t) = expected density of wolves from pack number 1
v(x, t) = expected density of wolves from pack number 2
p(x, t) = expected density of RLUs from wolf pack number 1
q(x, t) = expected density of RLUs from wolf pack number 2

1.1 Model for a Single Pack

In the simplest possible case, with only one wolf pack, we anticipate that two types of movements will dominate:

1. dispersal as the wolves search for food and other resources, and

2. movement back towards the den as the wolves return to the social organizing center (summer only).

Thus the word equation is:

Rate of change in expected wolf density

= Rate of change due to movement of wolves towards the den (summer only)
+ Rate of change due to dispersal of wolves away from high density

regions in search of food and other resources.

We model this situation by a nonlinear density-dependent diffusion-convection equation of the following form:

\frac{∂u}{∂t} = ∇ · [c_{u} \frac{(x − x_{u})}{|x − x_{u}|}u] + ∇ · [d_{u}(u)∇u] .  (1)

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 x_{u}. A related equation was proposed by Okubo (1986) to describe the swarming of insects (c.f. Murray, 1989).

1.2 Analytical Solutions for the Single-Pack Model

The simple single-wolf pack model (??) 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

d_{u}(u) = d_{0}u^{α}.  (2)

Integrating (??)–(??) twice using the conditions u → 0, ∂u/∂x → 0 as |x| → ∞ yields the steady state solution

u(x) = \{ \begin{matrix}(\frac{c_{u}}{d_{0}} )^{1/α} [\frac{U_{0}(1+α)}{2} − α|x − x_{u}|]^{1/α} , & |x − x_{u}| ≤ \frac{U_{0}(1+α)}{2α} \\ 0, & |x − x_{u}| > \frac{U_{0}(1+α)}{2α} \end{matrix}  (3)

(α > 0) where U_{0} is the total number of wolves

U_{0} = ∫_{−∞}^{∞} u(x) dx

(Figure 1) (c.f. Murray, 1989). When α = 0 the solution is

u(x) = \{ \begin{matrix}(\frac{U_{0}c_{u}}{d_{0}} )^{1/2} − \frac{c_{u}|x−x_{u}|}{d_{0}}, & |x − x_{u}| ≤ (\frac{U_{0}d_{0}}{c_{u}} )^{1/2} \\ 0, & |x − x_{u}| > (\frac{U_{0}d_{0}}{c_{u}} )^{1/2} . \end{matrix}  (4)

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 α = 0 (??)– (??) 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

u(x, t) = ∫_{−∞}^{∞} \frac{u(y, 0)}{2(πd_{0}t)^{1/2}} exp (\frac{−(x − x_{u} − y)^{2}}{4d_{0}t} ) dy.  (5)

Thus, for example, if all the wolves are initially at the den, then

u(x, 0) = U_{0}δ(x − x_{u})  (6)

and (??) describes a radially expanding, normally distributed population.

When c_{u} = 0 and α > 0 (??)–(??) becomes the porous media equation. Weak solutions for the radially symmetric case with the initial condition (??) can be calculated exactly by using a similarity transformation (see, for example, Murray (1989) or Grindrod (1991)). In the one-dimensional case, the solution,

u(x, t) = \{ \begin{matrix}\frac{1}{λ(t)} [1 − (\frac{x−x_{u}}{r_{0}λ(t)})^{2}]^{1/α} , & |x − x_{u}| ≤ r_{0}λ(t) \\ 0, & |x − x_{u}| > r_{0}λ(t), \end{matrix}  (7)

where

λ(t) = (\frac{t}{t_{0}} )^{1/(2+α)} , r_{0} = \frac{U_{0}Γ(\frac{1}{α} + \frac{3}{2})}{π^{1/2}Γ(\frac{1}{α} + 1)}, t_{0} = \frac{r_{0}^{2}α}{2d_{0}(α + 2)},

(Γ 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.

1.3 Wolf-Wolf Territorial Model

Inter-pack competition and subsequent territorial behavior is modelled by

1. changing equation (??) to include responses both familiar RLUs and to foreign RLUs from other packs, and

2. by including equations for the RLU densities themselves.

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:

Rate of change in expected density of wolves (pack 1)

= Rate of change due to movement of pack 1 wolves towards their den
+ Rate of change due to dispersal of pack 1 wolves
+ Rate of change due to movement of pack 1 wolves away from the RLUs

made by pack 2

+ Rate of change due to movement of pack 1 wolves towards the RLUs

made by pack 1.

\frac{∂u}{∂t} = ∇ · \{J_{c_{u}} + J_{d_{u}} + J_{a_{u}} + J_{b_{u}}\}  (8)

where the fluxes are given by

J_{c_{u}} = uc_{u}(q)\frac{(x−x_{u})}{|x−x_{u}|}, c_{u}(0) ≥ 0, \frac{dc_{u}}{dq} ≥ 0
J_{d_{u}} = d_{u}(u)∇u,  d_{u}(0) ≥ 0, \frac{dd_{u}}{du} ≥ 0
J_{a_{u}} = ua_{u}(q)∇q,  a_{u}(0) ≥ 0, \frac{da_{u}}{dq} ≥ 0
J_{b_{u}} = −ub_{u}(p)∇p, b_{u}(0) ≥ 0.

Here the rate of movement back to the den (which is situated at 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 dc_{u}/dq ≥ 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 x_{u} (c_{u}(0) > 0). The diffusive flux (d_{u}(u)) is similar to that described in the earlier simpler model (??).

Formulations for the fluxes away from competing RLUs (J_{a_{u}}) and towards familiar RLUs (J_{b_{u}}) assume that movement is in response to the gradient in expected RLU density (either ∇q or ∇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 J_{a_{u}} and J_{b_{u}} reflect that J_{a_{u}} describes movement away from competing RLUs and J_{b_{u}} describes movement 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:

Rate of change in expected density of wolves (pack 2)

= Rate of change due to movement of pack 2 wolves towards their den
+ Rate of change due to dispersal of pack 2 wolves
+ Rate of change due to movement of pack 2 wolves away from the RLUs

made by pack 1

+ Rate of change due to movement of pack 2 wolves towards the RLUs

made by pack 2.

\frac{∂u}{∂t} = ∇ · \{J_{c_{v}} + J_{d_{v}} + J_{a_{v}} + J_{b_{v}}\}  (9)

where

J_{c_{v}} = vc_{v}(p)\frac{(x−x_{v})}{|x−x_{v}|}, c_{v}(0) ≥ 0, \frac{dc_{v}}{dp} ≥ 0
J_{d_{v}} = d_{v}(v)∇v,  d_{v}(0) ≥ 0, \frac{dd_{v}}{dv} ≥ 0
J_{a_{v}} = va_{v}(p)∇p,  a_{v}(0) ≥ 0, \frac{da_{v}}{dp} ≥ 0
J_{b_{v}} = −vb_{v}(q)∇q, b_{v}(0) ≥ 0.

Lastly, we require equations to describe the changes in RLU density (p or q) with time. The equation for p is:

Rate of change in expected RLU density (pack 1)

= Low level of continual marking
+ Increased marking in the presence of RLUs from pack 2
+ First order decay of RLUs with time
\frac{∂p}{∂t} = u (l_{p} + m_{p}q) − f_{p}p,  (10)

and the equation for q is:

Rate of change in expected RLU density (pack 2)

= Low level of continual marking
+ Increased marking in the presence of RLUs from pack 1
+ First order decay of RLUs with time
\frac{∂q}{∂t} = v (l_{q} + m_{q}p) − f_{q}q.  (11)

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 Ω. A study of the system in two spatial dimensions requires Ω ⊂ IR^{2}. Later in this paper we will mathematically and numerically analyse the simpler one-dimensional case (Ω ⊂ IR).

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 Ω. This situation is described by zero-flux boundary conditions for u and v, namely

{J_{c_{u}} + J_{d_{u}} + J_{a_{u}} + J_{b_{u}}} · n on ∂Ω  (12)

and

{J_{c_{v}} + J_{d_{v}} + J_{a_{v}} + J_{b_{v}}} · n on ∂Ω,  (13)

where n is the outwardly oriented unit normal to the boundary of the solution domain, ∂Ω. Initial conditions, describing the expected spatial distributions of wolves and markings at the beginning of a study period, are given by

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

We now show that the zero-flux boundary conditions ((??)–(??)) guarantee a constant number of wolves for each pack within the domain Ω. At any given time, the total number of wolves from wolf pack 1 in the domain Ω is

∫_{Ω} u(x, t) dx.

The result follows directly from an application of the divergence theorem. For wolf pack 1

\frac{∂}{∂t} ∫_{Ω} u(x, t) dx = ∫_{Ω} \frac{∂u}{∂t}(x, t) dx

= ∫_{Ω} ∇ · {J_{c_{u}} + J_{d_{u}} + J_{a_{u}} + J_{b_{u}}} dx = ∫_{∂Ω} {J_{c_{u}} + J_{d_{u}} + J_{a_{u}} + J_{b_{u}}} · n ds = 0.

An analogous argument holds true for pack 2.

Noting that the area of Ω is given by

A = ∫_{Ω} dx,

we can calculate the average density of wolves from pack 1 throughout the region Ω as

U_{0} = \frac{1}{A} ∫_{Ω} u_{0}(x) dx  (15)

and the average density of wolves from pack 2 is

V_{0} = \frac{1}{A} ∫_{Ω} v_{0}(x) dx.  (16)

To simplify further analysis we nondimensionalize the model system (??)– (??) with boundary conditions (??)–(??) and intial conditions (??)). 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

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

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

b_{u}^{∗} = \frac{b_{u}V_{0}l_{q}}{L^{2}f_{p}^{2}} , b_{v}^{∗} = \frac{b_{v}V_{0}l_{q}}{L^{2}f_{p}^{2}} , μ_{p}^{∗} = \frac{m_{p}V_{0}l_{q}}{l_{p}f_{p}} , μ_{q}^{∗} = \frac{m_{q}U_{0}l_{p}}{l_{q}f_{p}} , ϕ = \frac{f_{q}}{f_{p}} .

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 Ω 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 (??)–(??) and

\frac{∂p}{∂t} = u (1 + μ_{p}q) − p.  (17)

\frac{∂q}{∂t} = v (1 + μ_{q}p) − ϕp.  (18)

The boundary conditions (??)–(??) remain unchanged, and an appropriate nondimensionalization of the initial data,

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

leaves the initial conditions (??) unchanged as well, once asterisks have been dropped. Observe that our nondimensionalization of space has rendered Ω an area (Ω ⊂ IR^{2}) or length (Ω ⊂ IR) equal to unity. Furthermore, with the nondimensionalization given here,

∫_{Ω} u(x, t) dx = ∫_{Ω} v(x, t) dx = 1  (19)

and thus, at any given time, u(x, t) and v(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.

1.4 Analytical and Numerical Results for the Territorial Model

We now analyze the two-wolf pack territorial model (??)–(??), (??)–(??), (??) and (??)–(??) mathematically and numerically. After making certain simplifying assumptions about the nonlinear functions in (??)–(??) 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 (??)–(??).

In particular we will consider the case where

\begin{matrix}c_{u}(·) = c_{v}(·) = c(·), & d_{u}(·) = d_{v}(·) = d(·), & a_{u}(·) = a_{v}(·) = a(·), \\ b_{u}(·) = b_{v}(·) = b(·), & μ_{p} = μ_{q} = μ, & ϕ = 1. \end{matrix}  (20)

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 onedimensional domain length equal to 1. The constraint (??) which was imposed by the nondimensionalization is thus re-written here in one spatial dimension as

∫_{0}^{1} u(x, t) dx = ∫_{0}^{1} v(x, t) dx = 1.  (21)

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

x_{u} = 0, x_{v} = 1.  (22)

Under these assumptions we will show that:

1. When avoidance of competing RLUs is modelled to include movement back towards the den or organizing center (J_{c_{u}} and J_{c_{v}} are nonzero and J_{a_{u}} and J_{a_{v}} are equal to zero (see (??)–(??))), 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.

2. 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 do not arise; solutions are spatially constant. (In this case, J_{c_{u}} and J_{c_{v}} are equal to zero and J_{a_{u}} and J_{a_{v}} are nonzero.) This result accentuates the importance of dens or organizing centers for the generation of territorial patterns.

1.4.1 Time-Independent Solutions With Dens or Organizing Centers

We now consider one-dimensional time-independent solutions to (??)–(??), (??)–(??), (??) and (??)–(??), (??)–(??), with

a(u) = a(v) = 0,  (23)

so that J_{a_{u}} and 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(·); that is, we choose

c(q) = cq, c(p) = cp (c constant).  (24)

Initially we consider the case where RLU marking is independent of the level of competing RLUs (μ = 0). We then relax this restriction and consider how patterns can evolve for nonzero μ.

\underline{The case μ = 0}

When μ = 0 the time-independent solutions to (??)–(??) are simply given by p = u and v = q. Thus steady state solutions to (??)–(??), (??)–(??) satisfy

0 = \frac{∂}{∂x} \{cuv + d(u)u_{x} − ub(u)u_{x}\}
0 = \frac{∂}{∂x} \{−cuv + d(v)v_{x} − vb(v)v_{x}\} .

Application of the zero-flux boundary conditions (??)–(??) yields

0 = cuv + d(u)u_{x} − ub(u)u_{x}  (25)
0 = −cuv + d(v)v_{x} − vb(v)v_{x}.  (26)

Note that these equations are symmetric under the transformation x → 1−x, u ↔ v. This symmetry, and the constraint (??), 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 (??)–(??) 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 ≥ 0. Defining

Γ(w) = ∫_{0}^{w} {d(ω) − ωb(ω)} dω,  (27)

we note that Γ(w) is a non-negative invertible function for w ≥ 0 because Γ(0) = 0 and Γ′(w) = d(w) − wb(w) > 0. Adding (??) and (??) and integrating with respect to x we observe that

Γ(u) + Γ(v) = Γ(u(0)) + Γ(v(0)) = k (constant)  (28)

and thus (??) is conserved for all u(x) and v(x) satisfying (??)–(??). We use the fact that Γ is invertible to separate u and v in (??)–(??) and rewrite them as

0 = cuΓ^{−1}(k − Γ(u)) + (d(u) − ub(u))u_{x}  (29)
0 = −cvΓ^{−1}(k − Γ(v)) + (d(v) − vb(v))v_{x}.  (30)

Integration of (??) yields u(x) and similarly integration of (??) yields v(x). By way of examples, we now consider two special cases:

1. 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 (??) yields Γ(w) = dw and thus Γ^{−1}(w) = d^{−1}w. Therefore equation (??) indicates that u + v = k/d is conserved across the solution domain 0 ≤ x ≤ 1. Equations (??) and (??) become

u_{x} = \frac{−c}{d}u (\frac{k}{d} − u)  (31)
v_{x} = \frac{c}{d}v (\frac{k}{d} − v) .  (32)

These are simply versions of the logistic equations with space as the independent variable. Explicit solutions for u and v are given by

u(x) = \frac{ku(0) exp\{−\frac{d^{2}}{ck}x\}}{k + du(0)[exp\{−\frac{d^{2}}{ck}x\} − 1]}  (33)

v(x) = \frac{kv(0) exp\{\frac{d^{2}}{ck}x\}}{k + dv(0)[exp\{\frac{d^{2}}{ck}x\} − 1]}  (34)

where k is k = d(u(0) + v(0)) (see (??)). By integrating (??) and (??) with respect to x and applying the constraint (??) we have

u(0) = (u(0) + v(0))\frac{exp \{−\frac{d}{c(u(0)+v(0))^{2}} \} − 1}{exp \{−\frac{d}{c(u(0)+v(0))}\} − 1}  (35)

v(0) = (u(0) + v(0))\frac{exp \{\frac{d}{c(u(0)+v(0))^{2}} \} − 1}{exp \{\frac{d}{c(u(0)+v(0))}\} − 1}.  (36)

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 (??)–(??) where k = d(u(0) + v(0)). A numerical solution to the related time-dependent problem ((??)–(??), (??)–(??), (??) and (??)–(??), (??), (??), (??), a(w) = 0, d(w) = d, b(w) = 0, and μ = 0) achieved the predicted heterogeneous steadystate pattern for large time (Figure 2). Initial conditions were given by

u(x, 0) = 1 − H(1/2), v(x, 0) = H(1/2).  (37)

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

(a) The total expected wolf density (u+v) remains constant throughout and is not reduced near the inter-pack boundary and

(b) 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 μ ≠ 0.

Before considering the μ ≠ 0 case, we briefly analyse a case with μ = 0 and nonlinear diffusion.

2. d(w) = dw (d constant), and b(w) = b (b constant), d > b:

Here equation (??) yields Γ(w) = (d − b)w^{2}/2 and thus Γ^{−1}(w) = \sqrt{2w/(d − b)}. Therefore equation (??) says that u^{2} + v^{2} = 2k/(d − b) is conserved across the solution domain 0 ≤ x ≤ 1. Equations (??) and (??) become

0 = u [c (\frac{2k}{d − b} − u^{2})^{1/2} + (d − b)u_{x}]  (38)

0 = v [−c (\frac{2k}{d − b} − v^{2})^{1/2} + (d − b)v_{x}] .  (39)

Integration of (??)–(??) and application of the integral constraint (??) yield the following solution which is valid when d − b ≥ 2c/π:

u(x) = \{ \begin{matrix}\frac{2}{1+(d−b)(2−π/2)/c} & 0 ≤ x ≤ x_{0} \\ \sqrt{\frac{2k}{d−b}} cos (\frac{c(x−x_{0})}{d−b}) & x_{0} < x < x_{1} \\ 0 & x_{1} ≤ x ≤ 1 \end{matrix}  (40)

v(x) = \{ \begin{matrix}0 & 0 ≤ x ≤ x_{0} \\ \sqrt{\frac{2k}{d−b}} sin (\frac{c(x−x_{0})}{d−b}) & x_{0} < x < x_{1} \\ \frac{2}{1+(d−b)(2−π/2)/c} & x_{1} ≤ x ≤ 1 \end{matrix} ,  (41)

where

x_{0} = \frac{1}{2} (1 − \frac{(d − b)π}{2c} ) , x_{1} = \frac{1}{2} (1 + \frac{(d − b)π}{2c} ) .

The solution for the case 0 < d − b < 2c/π is given by

u(x) = \sqrt{\frac{2k}{d − b}} cos (\frac{c(x − x_{0})}{d − b})

v(x) = \sqrt{\frac{2k}{d − b}} sin (\frac{c(x − x_{0})}{d − b}.)

This solution is shown in Figure 3. We observe that the total wolf density is, in fact, higher near the pack boundary than it is near the den.

While such results clearly are not commensurate with field observations of wolves discussed elsewhere, we have included the solution primarily because it provides insight into the effect of nonlinear diffusion on the formation of territorial patterns. Observable effects are:

(a) The formation of abrupt boundaries (with discontinuous u_{x} and v_{x}); the steady state solution is thus a weak solution to the original time dependent problem (see the earlier section for a similar result for the one-equation model).

(b) In this case, an increase in the overlap of adjacent wolf territories.

\underline{The case μ ≠ 0}

When μ ≠ 0 the time-independent solutions to (??)–(??) are given by

p = \frac{u(1 + μv)}{1 − μ^{2}uv}  (42)
q = \frac{v(1 + μu)}{1 − μ^{2}uv}  (43)

with spatial derivatives given by

p_{x} = \frac{1}{(1 − μ^{2}uv)^{2}} [u_{x}(1 + μv) + v_{x}μu(1 + μu)]  (44)
p_{x} = \frac{1}{(1 − μ^{2}uv)^{2}} [v_{x}(1 + μv) + u_{x}μv(1 + μu)] ,  (45)

respectively. The possibility of unbounded p and q thus emerges for the case 1 − μ^{2}uv = 0. We therefore restrict our discussion to cases where μ is sufficiently small so that 1 − μ^{2}uv > 0. Thus steady state solutions to (??)–(??), (??)–(??) satisfy

0 = \frac{∂}{∂x} \{cuq + d(u)u_{x} − ub(p)p_{x}\}
0 = \frac{∂}{∂x} \{−cvp + d(v)v_{x} − vb(q)q_{x}\}

with p, q, p_{x} and q_{x} defined by (??)–(??). Application of the zero-flux boundary conditions (??)–(??) yields

0 = cuq + d(u)u_{x} − ub(p)p_{x}  (46)
0 = −cup + d(v)v_{x} − vb(q)q_{x}  (47)

Thus (??)–(??) defines a system of two coupled nonlinear ordinary differential equations for u and v. Therefore, given any prescribed d(·) and b(·) 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(·) and b(·). Substitution of (??) and (??) into (??)–(??) yields

0 = c\frac{uv(1 + μu)}{1 − μ^{2}uv} + du_{x}  (48)
0 = −c\frac{uv(1 + μv)}{1 − μ^{2}uv} + dv_{x}  (49)

We can now formulate a conservation law as we did earlier for the μ = 0 case (see (??)). Defining

Γ_{1}(w) = ∫_{0}^{w} \frac{dω}{1 + μω} = log(1 + μw)/μ

we observe from (??)–(??) that

Γ_{1}(u) + Γ_{1}(v) = Γ_{1}(u(0)) + Γ_{1}(v(0)) = k (constant).

Thus

(1 + μu)(1 + μv) = exp(μk)  (50)

describes u in terms of v and vica-versa.

Substitution of (??) into (??) yields

0 = −c\frac{v (\frac{exp(μk)}{1+μv} − 1) /μ}{1 − μv (\frac{exp(μk)}{1+μv} − 1)} + \frac{dv_{x}}{1 + μv}  (51)

while substitution into (??) yields a similar equation for u. Equation (??) is separable, and integration yields an implicit formula for v(x), namely

\frac{v(x)}{(E − μv(x))(1 + μv(x))^{E}} = \frac{v(0) exp\{\frac{cEx}{dμ}\}}{(E − μv(0))(1 + μv(0))^{E}} ,

where we have written E = exp(μk) − 1.

For the special case E = 1 we can calculate the solution explicitly as

v(x) = \frac{1}{2μ^{2}} [(μ^{2}v(0)^{2} − 1) exp \{\frac{−cx}{dμ}\} /v(0)

+ ((μ^{2}v(0)^{2} − 1)^{2} exp \{\frac{−2cx}{dμ}\} /v(0)^{2} + 4μ^{2})^{1/2}] . (52)

Noting the symmetry of (??)–(??) under the transformation x → 1 − x and u ↔ v we can deduce that u(0) = v(1). From equation (??) we observe that, when E = 1, (1 + μu(0))(1 + μv(0)) = 2, and thus

(1 + μv(1))(1 + μv(0)) = 2.  (53)

Integration of (??) with respect to x and application of the constraint (??) defines v(0) in terms of c/d, and μ. Specifically, v(0) must satisfy

1 = −\frac{dv(0)}{2μc} [1 − exp \{\frac{−c}{dμ}\}] + \frac{d}{c} [− ((\frac{v(0)}{2μ})^{2} exp \{\frac{−2c}{dμ}\} + 1)^{1/2}

−\frac{1}{2} log (1 − ((\frac{v(0)}{2μ})^{2} exp \{\frac{−2c}{dμ}\} + 1)^{−1/2})

+\frac{1}{2} log (1 + ((\frac{v(0)}{2μ})^{2} exp \{\frac{−2c}{dμ}\} + 1)^{−1/2})

+ ((\frac{v(0)}{2μ})^{2} + 1)^{1/2} + \frac{1}{2} log (1 − ((\frac{v(0)}{2μ})^{2} + 1)^{−1/2})

−\frac{1}{2} log (1 + ((\frac{v(0)}{2μ})^{2} + 1)^{−1/2})] .  (54)

We have thus derived a one-parameter family of solutions to (??)–(??) for the case E = 1; when given a value for μ, we can calculate v(0) and c/d by simultaneously solving two equations, the first of which is equation (??) in conjunction with (??) for the definition of v(1) in terms of v(0), and the second of which is equation (??).

In Table 1 we computed the values of v(0) and c/d, while varying μ. The corresponding profiles for v(x), as given by (??), 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 μ = 0 case.

A numerical solution of the related time-dependent problem ((??)–(??), (??)–(??), (??) and (??)–(??), (??), (??), (??), a(w) = 0 and b(w) = 0) for the μ = 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 (??). (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 do not have such a maximum (Figure 6). We can show that this must indeed be the case by looking at equations (??) and (??) in conjuction with equation (??). By differentiating (??) with respect to x we observe that

0 = \frac{v_{x}}{1 + μv} + \frac{u_{x}}{1 + μu}

and thus (??) and (??) can be simplified to give

p_{x} = \frac{(1 + μu)(μu − 1)}{(1 − μuv)^{2}}v_{x}  (55)
q_{x} = \frac{(1 + μv)(μv − 1)}{(1 − μuv)^{2}}u_{x}.  (56)

Thus we can calculate the critical points for p and q as

p_{x} = 0 ⇔ v_{x} = 0 or u = 1/μ  (57)
q_{x} = 0 ⇔ u_{x} = 0 or v = 1/μ.  (58)

Equations (??) and (??) indicate that the signs of u_{x} and v_{x} do not change. (You may recall our earlier assumption that 1 − μ^{2}uv > 0). Thus an interior maximum for p requires that u = 1/μ and an interior maximum for q requires that v = 1/μ. Substitution of u = 1/μ and E = exp(μk) − 1 into (??) indicates that an interior maximum for p can only occur if

u = \frac{1}{u}, v = \frac{E − 1}{2μ}  (59)

and substitution of v = 1/μ and E = exp(μk) − 1 into (??) indicates that an interior maximum for q can only occur if

u = \frac{E − 1}{2μ}, v = \frac{1}{u}.  (60)

Thus for the case E = 1, an interior maximum for p requires that v = 0, however never occurs (see equation (??)). 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 (??) and (??), 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 (??) evaluated at (??), namely

pmax = \frac{1 + E}{μ(3 − E)}.  (61)

Similarly an interior maximum for q would be given by (??) evaluated at (??), namely qmax = pmax. If bounded interior maximums p and q were to occur, they would thus be for E in the range 1 < E < 3.

Our approach to investigating this possibility was to solve the full nonlinear time-dependent problem numerically and look for a steady state solution for large time. Numerical solutions did achieve a spatially heterogeneous steady state profile for large time. A sample result for c = 1, d = 0.3, μ = 1.1 is shown in Figure 7. Note that distinct maximums do occur for p and q. The numerically calculated value for E was E = 2.3 and thus the predicted values for pmax and qmax were 4.3. These compare well with the computed results. According to (??) and (??) the maximum for p should occur at u = 1/μ = 0.91 and v = (E−1)/(2μ) = 0.59 and the maximum for q should occur at v = 1/μ = 0.91 and u = (E − 1)/(2μ) = 0.59. These predictions are also borne out by the numerical solution (Figure 7). This numerical solution encapsulates two key biological observations about territorial pattern formation in wolves (discussed elsewhere), namely the decreased wolf density and the increased RLU densities near territorial boundaries.

1.5 Time-Independent Solutions With no Den or Or-

ganizing Center

Jim, maybe we should talk about this when I visit.

List of Figures

Figure 1: One-dimensional solutions to the nonlinear density-dependent diffusionconvection equation (??)–(??) are shown schematically. (a) A heterogeneous steady-state solution is shown for α = 1/2 (see (??). (b) A radially expanding solution is shown for c_{u} = 0 (porous media equation) (see (??)).

Figure 2: Numerical solution to the time-dependent problem achieves a steady-state solution for large time. Shown is the solution to the timedependent problem: (??)–(??), (??)–(??), (??) and (??)–(??), (??), (??), (??), a(w) = 0, d(w) = d, b(w) = 0, and μ = 0 for t =??. Initial conditions are given by (??). The Method of Lines and Gear’s Method were used to solve the system.

Figure 3: The analytical solution to (??)–(??) is given by (??). The solution is shown here for b =?, c =? and d =?.

Figure 4: The analytical solution to (??), (??) for E = 1 is given by (??). Solutions are shown here for the values given in Table 1, with μ varying from 0.42 (shallow curve) to 0.48 (steep curve).

Figure 5: Cumulative wolf-pack densities u(x) + v(x) and cumulative RLU densities p(x) + q(x) were calculated using (??), u = v(1 − x), (??) and q = p(1 − x). The values for μ is μ = 0.48 and the values for v(0) and c/d are taken from Table 1. Note that the cumulative pack density decreases and cumulative RLU density increases near x = 1/2.

Figure 6: A numerical solution the time-dependent problem achieves a steadystate solution for large time Shown is the solution to the time-dependent problem: (??)–(??), (??)–(??), (??) and (??)–(??), (??), (??), (??), a(w) = 0 and b(w) = 0 for the μ = 0.48 case from Table 1 (c = 1 and d = 1/7.3584) for t =??. Initial conditions were given by (??). The Method of Lines and Gear’s Method were used to solve the system.

Figure 7: A numerical solution the time-dependent problem achieves a steadystate solution for large time Shown is the solution to the time-dependent problem: (??)–(??), (??)–(??), (??) and (??)–(??), (??), (??), (??), a(w) = 0 and b(w) = 0 for the case μ = 1.1, c = 1, d = 0.3. Initial conditions were given by (??). (a) Profiles for u(x), v(x), p(x) and q(x). Note that both p(x) and q(x) have interior maximums near the point x = 1/2. (b) Profiles for u(x), v(x), the cumulative wolf pack density (u(x) + v(x)), and E + 1 = exp(μk) = (1 + μu(x))(1 + μv(x)). Note that the cumulative wolf pack density is lowest at x = 1/2, and E + 1 is constant across the domain. The Method of Lines and Gear’s Method were used to solve the system.