The domestication of prehistoric wolves

Alice Henman teaches us to sit…, stay…, wait…wait…no, come back! Fenton!

post

Known as man’s best friend, dogs have evolved alongside humans for thousands of years and come in many shapes and forms, aiding mankind in a huge range of specialised roles. There are herding and guarding breeds to help manage and protect livestock; hunting dogs used for tracking, killing or retrieving game; dogs used in the workforce for their powerful noses or strong builds; and service dogs for assisting people with disabilities. For centuries, we even had a breed called the turnspit, specialised for turning a wheel, often to power kitchen appliances. But how did we get here?

Turnspit dogs ran in wheels to turn meat, depicted here in a Welsh pub in 1797.

It is well known that we’ve been selectively breeding dogs for these roles for centuries, but how we came into possession of such docile animals in the first place is still up for debate. We know that dogs evolved from prehistoric grey wolves, while grey wolves themselves are still around today. At some point, an evolutionary branch occurred, separating the wolves that remained wild from those that would eventually become our household pets. Although it is not known exactly what caused this distinction, there have been many hypotheses.

The domestication process. (Source)

Some theories suggest that complementary hunting styles led to collaboration between our respective prehistoric species, and the relationship grew from there. Others believe that humans may have hand-reared wolf pups from an early age; the pups would imprint on their carer, forming the initial bond. Aggressive wolves would be culled, whereas successfully socialised wolves would be kept and bred to produce increasingly domesticated companions. The final popular theory is of proto-domestication: the idea that early wolves domesticated themselves in favour of scavenging human scraps, which would’ve been more reliable as a food source than wild prey.

We can model the evolution of these prehistoric wolves into our modern-day pets and their cousins, the grey wolf, using mathematics! In this case we will test the mathematical validity of the proto-domestication theory using a method called adaptive dynamics.

Adaptive dynamics

Adaptive dynamics takes the theory of natural selection and applies it to the parameters defining the dynamical systems that govern a population. This means we assume that evolution will favour traits that are advantageous to the population. In this case, we assume a mutant trait will be successful if it has a faster ‘growth rate’ relative to the trait already resident in the population.

Before we go any further, there are some key simplifications at play in adaptive dynamics that must be stated.

  1. Mutations are small: a mutant individual doesn’t differ much from its parent.
  2. Mutations are rare: The number of mutants that arise at any one time is very small. This allows us to approximate the total population with the resident population (the one in place before any mutations).
  3. The evolutionary timescale is much larger than the dynamical system timescale. The population is allowed to reach equilibrium between successive mutations.

Generally, adaptive dynamics considers the evolution of a single evolutionary trait. Backed up by Dmitri Belyaev’s studies on the domestication of foxes, begun in the 1950s, we can reasonably take the trait controlling the domestication of wolves to be friendliness, tameness, or something similar. The study also found that many traits found in domesticated animals, like ear floppiness, softened features, shorter snouts and smaller bodies, were linked, and often emerged on their own among more docile individuals in the species. This makes the transition from wolves to domestic dogs easier to model, as some characteristics don’t need to be accounted for individually, and instead we can just focus on the single friendliness trait.

Beware of the $f=0$

We can denote the friendliness attribute by the parameter $f\in[0,1]$ in the dynamical system we use to describe the prehistoric wolf population. A friendliness value of $0$ corresponds to completely wild wolves with a high aversion to humans. They are therefore strong hunters who don’t benefit from scavenging leftovers. An $f$ value of 1, however, marks a fully domesticated wolf relative to the standard at the time, before selective breeding. These wolves will be the least fearful of humans and are therefore able to benefit from scavenging scraps from human settlements.

Since we are dealing with wolf populations, it makes sense to use a classic Lotka–Volterra model, often used for predator–prey type relationships, but instead the focus will be on interspecies interactions rather than between wolves and their prey.

For two sub-populations of wolves with different friendliness, we say that the first sub-population, $W$, grows proportionally to its current size, but is limited by competition from itself and the second sub-population, $W’$:
\begin{equation*} \frac{\mathrm{d}{W}}{\mathrm{d}t} = \underbrace{rW}_{{\text{growth}}} – \underbrace{\frac{r}{K}W^2}_{{\substack{\text{self-}\\\text{competition}}}} – \underbrace{cWW’}_{{\text{with }W’}}. \end{equation*}
Here, $r$ is the usual growth rate, encoding information about breeding success and litter sizes; and $K$ is the carrying capacity of the sub-population of wolves $W$, ie their limiting population size in the absence of any other sub-populations.

More generally, for many sub-populations $W(f,t)$ with different values of $f$ this extends to (after some relabelling)
\begin{equation*} \frac{\mathrm{d}W(f,t)}{\mathrm{d}t} = rW(f,t)\left[1\!-\!\sum_{f’}\frac{\alpha(f-f’)W(f’,t)}{K(f)}\right]. \label{pop_density} \end{equation*}
The term $\alpha(f-f’)$ is a Gaussian function describing the competition a wolf of friendliness $f’$ exerts on another with trait value $f$,
\begin{equation*} \alpha(f-f’) = \exp{\left(\frac{-(f-f’)^2}{2\sigma_\alpha^2}\right)}. \label{branch_comp} \end{equation*}
This competition distribution corresponds to individuals with similar traits competing more aggressively, while individuals with sufficiently different trait values compete less strongly, where the spread of the competition is characterised by the width, $\sigma_\alpha$. This is because wolves with the same or similar tameness use similar resources: wild wolves compete for the prey they are hunting, while docile ones compete for human scraps. Two examples can be seen in the graph below.

We can use Gaussian distributions to represent the competition exerted by the trait $f’$ on the trait $f$ by the difference $f-f’$. The two widths, $\sigma_\alpha$, represent different spreads of trait values.

For the carrying capacity $K(f)$, we choose a double Gaussian to capture the potential benefits of both extremes. As seen below, the two resource components are centred at $f=0$ (wild) and $f=1$ (docile) and are described by the function
\begin{equation*} K(f) = \underbrace{A\exp{\left(\frac{-f^2}{2\sigma_\text{w}^2}\right)}}_{{\text{wild resource}}} + \underbrace{B\exp{\left(\frac{-(f-1)^2}{2\sigma_\text{d}^2}\right)}}_{{\text{human resource}}}. \end{equation*}

The carrying capacity of the wolf population varies with trait value $f$ and the niche widths $\sigma_\text{w}$ and $\sigma_\text{d}$. Each hump indicates trait values that promote the highest equilibrium population densities. Shading excludes impossible values of $f$. Here, $A=2.5$ and $B=2$.

When these widths are small enough, the Gaussians representing each food source are more distinct and two separate peaks are visible; in other words, it’s better to be an expert than a jack of all trades. The wildest wolves will be the best at hunting prey: as aggression diminishes, this ability is lessened, so friendlier wolves will have less hunting success. Conversely, the most fearless wolves will get first dibs on human scraps, leaving the more cautious wolves with less.

Similarly, the idea of wild wolves being the best hunters can be flipped: those who have less success with wild prey are more likely to turn to self-domestication.

Larger spreads produce a significant overlap between the two resources, leading to the loss of these peaks and the emergence of some optimum middle ground between wild and human resources.

Who’s a good wolf?

Now that we have our model set up, we can get started on the adaptive dynamics analysis. First off, the resident trait value will always be denoted $f=x$; this is the value of $f$ before it mutates at any given point, while $f=y$ will be used as the placeholder for each mutant friendliness value.

As per our assumptions of adaptive dynamics, mutant rarity tells us that initially we can neglect the mutant population density, ie $W’ = W(y,t)\approx0$, and a large evolutionary timescale lets us set $W(x,t)$ to its equilibrium value, $W^*(x)$. The equilibrium occurs where $\frac{\mathrm{d}W(x,t)}{\mathrm{d}t}=0$; since we assume the population is in equilibrium before mutations arise, our equation for $\frac{\mathrm{d}W}{\mathrm{d}t}$ says this is equivalent to
\begin{align*} rW^*(x)\left[1-\frac{\alpha(x-x)W^*(x)}{K(x)}\right] &= 0 \\ \implies W^*(x) &= K(x). \label{eqn:residentdensity} \end{align*}
The first thing we can do with this model is test the success of a mutant wolf (one which is more or less human resistant than its pack) when the rest of the population all shares the same friendliness value. We do this by finding what The hitchhiker’s guide to adaptive dynamics calls the invasion fitness: the expected per capita growth rate of a rare mutant in the environment set by the resident population.

Denoted $F_x(y)$, where $x$ is the resident trait and $y$ is the mutant trait, this is
\begin{align*} F_x(y) &= \left.\frac{1}{W(y,t)}\frac{\mathrm{d}W(y,t)}{\mathrm{d}t}\right|_{\begin{subarray}{l}W(x,t)=W^*(x)=K(x),\\W(y,t)=0\end{subarray}}\\ &= r\left[1-\frac{\alpha(y-x)K(x)+\alpha(x-x)\cdot 0}{K(y)}\right]\\ &= r\left[1-\frac{\alpha(y-x)K(x)}{K(y)}\right]. \end{align*}
If the invasion fitness of $y$ is positive, the mutant can invade; if it is negative, it will die out.

Given our carrying capacity $K(x)$ above, the plot below shows $F_x(y)$ for $x=0.2$ (pink), $0.5$ (purple) and $0.8$ (blue): resident traits of mostly wild, intermediate and mostly docile.

Notice that $F_x(x)=0$ always: a mutant identical to the resident is neither advantaged nor disadvantaged. So, given the resident trait $x$, which direction for $y$ is evolutionarily more successful? The one that’s positive!

Given the carrying capacity, $K(f)$, shown above with $\sigma_\text{w}=1,\sigma_\text{d}=0.5$ (orange), here are three plots of the invasion fitness, $F_x(y)$, for different resident traits, $x$. Notice that $F_x(x)=0$ in every case: if the mutant is the same as the resident trait, there will be no (dis)advantage. The strength of competition width is $\sigma_\alpha=0.5$.

So the sign of the invasion fitness, $F_x(y)$, matters. This can also be visualised on something called a pairwise invasibility plot (PIP) which shows the sign of $F_x(y)$ on a plot of resident trait value by mutant trait value. These are plotted below.

Plots of mutant trait value by resident trait value. Green regions show where the invasion fitness is positive; white regions are negative. The regions are separated by black contours of $F_x(y)=0$. The left-hand plot shows that the curves of fixed $x$ from the previous graph exhibit the same invasion fitness polarity. The right-hand plot shows how traits will evolve along the diagonal: if a mutant trait is in the green region it will be successful and the population will evolve to adopt this new trait. But if it lies in the white region, the resident trait will remain the same. This evolution sequence is shown for a limited number of mutations but the pattern shows evolution will continue towards the point where the contours intersect. Parameters: $A=2.5$, $B=2$, $\sigma_\text{w} = 1$, $\sigma_\text{d} = 0.5$ and $\sigma_\alpha = 0.6$.

The selection gradient, $s(x)$, for the resident trait, $x$, tells us whether friendlier or less friendly mutants are favoured for each $x$. This plot is based
on the carrying capacity from above with $\sigma_\text{w} = 1$, $\sigma_\text{d} = 0.5$.

We have assumed that mutations are small, so another way to find which direction for $y$ is more successful is to look at the slope of the graph at $y=x$, known as the selection gradient:
\begin{align*} s(x) &= \left.\frac{\partial F_x(y)}{\partial y}\right|_{y=x} \\ &=r\,\Biggl[-\underbrace{\left.\frac{\partial \alpha(y-x)}{\partial y}\right|_{y=x}}_{=\;0} + \frac{\partial K(x)}{\partial x}\frac{1}{K(x)}\Biggr]\\ &=\frac{r}{K(x)}\frac{\partial K(x)}{\partial x}. \end{align*}
The selection gradient tells us information about the speed and direction of evolution: positive $s$ means slightly friendlier mutants are favoured; negative $s$ means slightly less-friendly ones are favoured.

This also tells us whether we move up or down the diagonal of the PIP: if it is green above the diagonal and white below, a friendlier wolf will be more successful so $s(x)>0$ (we move upward and rightward); but the converse means $s(x) < 0$ and we move downward and leftward. To the left you can see $s(x)$ for different values of the resident trait, $x$.

A golden retri-evolutionarily singular strategy

The trait value at which evolution gets interesting is where $s(x^*)=0$: neither friendly nor unfriendly is favoured, so this trait value is called the evolutionarily singular strategy. This can be seen visually as the point at which the $F_x(y)=0$ contours cross in the PIP.

‘Gets interesting’ can mean many things, but in our case, the point of interest we are looking for is where the singular strategy is additionally a branch point. Branch points are defined by their ability to be both convergent and unstable evolutionary waymarks: evolution is attracted to the point but selection can then drive the population to split into distinct branches which diverge away from one another towards different values of the trait, $f$.

Convergence comes from requiring $s'(x^*) < 0$; instability thereafter requires that $F_{x^*}(y)$ has a local minimum at $y=x^*$ so that any mutation away from $x^*$ grows. Some mathematical manipulation can show that meeting these conditions corresponds to satisfying the inequality \begin{alignat*}{3} \left.\frac{\partial^2 F_x(y)}{\partial x^2}\right|_{y=x=x^*} &>\;& \left.\frac{\partial^2 F_x(y)}{\partial y^2}\right|_{y=x=x^*} &>0. \end{alignat*} In our case, this implies
\begin{equation*}\begin{split} r\left(\frac{1}{\sigma_\alpha^2} – \frac{1}{K(x^*)}\left.\frac{\partial^2 K(x)}{\partial x^2}\right|_{x=x^*}\right) \hspace{2cm} \\ > r\left(\frac{1}{\sigma_\alpha^2} + \frac{1}{K(x^*)}\left.\frac{\partial^2 K(x)}{\partial x^2}\right|_{x=x^*}\right) >0. \end{split} \end{equation*}
Let’s test it out on the wolf population. Firstly, our evolutionary singular strategies can be found simply from solving $\frac{\partial K(x)}{\partial x}=0$ (the stationary points of the carrying capacity), which can be done numerically.

Benji’s favourite branch point is under one, in the shade. (Benji, $f=1$, will happily eat human scraps.)

Next, it’s clear from the inequality above that we must have
\begin{equation*} \left.\frac{\partial^2 K(x)}{\partial x^2}\right|_{x=x^*} < 0, \end{equation*} to guarantee convergence, so branch points must be maxima of $K(f)$. The instability condition is then only satisfied for \begin{align*} \frac{K(x^*)}{\sigma_\alpha^2} > &-\left. \frac{\partial^2 K(x)}{\partial x^2}\right|_{x=x^*} = \left| \frac{\partial^2 K(x)}{\partial x^2}\right|_{x=x^*}\\ \implies&\frac{K(x^*)}{\hspace{2em}\left|\frac{\partial^2 K(x)}{\partial x^2}\right|_{x=x^*}} > \sigma_\alpha^2. \end{align*}

Since both sides are positive, this can hold true if the wolf–wolf competition has a narrow enough width! There are two possible scenarios where this can happen — look again at the plot of $K(f)$ above.

Case 1: The yellow plot of $K(f)$:

  • The niche widths are narrow.
  • The branch point resides near $f=0$.

This is the branch of interest to our domestication story because initially all wolves are wild and use the wild resource.

Case 2: The orange plot of $K(f)$:

  • The niche widths are wide.
  • The branch point is at a mid-friendliness value.

Both cases can be shown in stochastic simulations under the appropriate constraints derived for branching earlier. You can see these on the next page.

From the plots, it seems as though case 1 is the more realistic choice, although it is worth noting that our example of case 2 is an extreme one. By altering the values of $A$, $B$, $\sigma_\text{w}$ and $\sigma_\text{d}$, the branch point can be pushed closer to $0$ to be more representative of where we would expect the true branch to lie.

In fact, we would expect to find $A$ much larger than $B$ in the early beginnings of wolf proto-domestication (corresponding to the purely wild wolf population having a much higher capacity than the purely ‘domestic’) and when this is true, the branch point is found to be very near zero despite the case we are in.

This is good news since grey wolves are essentially the direct descendants of pre-historic wolves, so it makes sense that the branch point lies right on $f=0$ and dogs have diverged from wolves at that point leaving the wolves unperturbed, rather than all wolves evolving as a whole towards a more docile state before going back on themselves (evolutionarily speaking) after branching.

Stochastic simulations of friendliness evolution in prehistoric wolves over time. Values of $\sigma_\alpha$ have been chosen to satisfy the conditions for a branch point near $0$ in both cases. Each coloured curve represents a different realisation of the simulation to show the pattern is not unique. In both cases, one branch would evolve to consist of purely wild wolves, while the other would grow in friendliness until they eventually become domestic enough for selective breeding to take over (which would require a whole other model). Parameters: $A=2.5$, $B=2$, maximum mutational step $=0.01$. Case 1: $\sigma_\text{w} = 0.5$, $\sigma_\text{d}=0.25$, $\sigma_\alpha = 0.4$. The simulation shows that in the presence of humans, the population would branch at some value of $f$ very close to zero. Case 2: $\sigma_\text{w} = 1$, $\sigma_\text{d}=0.5$, $\sigma_\alpha = 0.6$. The simulation shows that in the presence of humans, the population would grow in friendliness until some critical value is reached, at which point the population splits into increasingly wilder and friendlier wolves. Images: Grey wolf image, Seney Natural History Association, CC BY-SA 2.0.

Give paw

So what have we found out? For the model specified at the beginning of this article, the maths supports the proto-domestication theory that the benefit of human scraps could have driven the evolution of prehistoric wolves into dogs and grey wolves.

What’s most fascinating is not just that self-domestication was attractive enough to cause this change, but that the wild wolf population remained. Evolution was not unidirectional; instead, it supported the creation of two distinct populations within the same geographical area! This is important as speciation is often considered to be driven by geographical separation which consequently causes a species to evolve to adapt to two different environments, which wouldn’t have necessarily been the case for these wolves. Furthermore, this is an important result as it shows the usefulness of mathematics in situations where we can’t know the answer to a question using data or experiments (in this case because it happened so long ago).

Although we still don’t have a concrete answer for the origin of our four-legged friends this maths acts as one more piece of the puzzle in the quest to find out.

Alice is starting a PhD in Durham in mathematical biology (on metals and microbes). She is part of her university college’s grow society and ‘less culty’ volleyball society. Just don’t ask her to look after your plants or actually win a volleyball game.

More from Chalkdust