editorial note: These notes pre-date the formal start of my online laboratory notebook, Feb 2 2010: The Lab Notebook Goes Open and were adapted from a LaTeX document in which I kept notes on this topic during my summer at IIASA. Lacking a proper notebook then, documents like this one were updated periodically and occassionally branched into new ones. The post date represents the last time the LaTeX document was edited in the course of that research.
We consider the classic adaptive dynamics model for evolutionary branching under competition for a limiting resource by Dieckmann and Doebeli (1999). The population dynamics of \(N(x,t)\) individuals each with focal trait \(x\) are given by a continuous time birth-death process with rates
\[\begin{aligned} b(N) = & r \nonumber \\ d(N) = & r \frac{\sum_yN(y,t)C(x,y)}{K(x)} .\label{general_logistic}\end{aligned}\]
The focal trait \(x\) can be thought of as some continuous phenotypic variable, such as beak size, which affects resource consumption. Niche differences driving speciation arise from differences in this trait. In the model, \(K(x)\) is the resource kernel which determines the equilibrium population density and \(C(x,y)\) is a the competition kernel, describing the relative change in death rate of individuals of type \(x\) due to competition by individuals of type \(y\). Given \(C(x,x)=1\), the equilibrium density of a monomorphic population is \(N^{\ast}(x)=K(x)\). We assume the following Gaussian forms,
\[K(x)=K_0e^{-x^2/(2\sigma_k^2)}, \label{K}\]
and
\[C(x,y)=e^{-(x-y)^2/(2\sigma_c^2)}, \label{C}\]
where \(\sigma_k\) and \(\sigma_c\) are scale factors for the resource distribution and competition kernels respectively. \(\sigma_c\) can be though of as the niche size. We focus on the case \(\sigma_k>\sigma_c\), for which the model has an evolutionary branching point at \(x=0\)
The evolutionary dynamics are specified only at the phenotypic level: We assume that with probability \(\mu\) an individual birth results in a mutant offspring, and that the mutant trait is chosen from a Gaussian distribution centered at the trait value of the parents and with a variance \(\sigma_{\mu}^2\).
Classically, the approach to the branching point is governed by the canonical equation of adaptive dynamics. To account for branching we have to revisit the original derivation of this equation where we relax the limit of small mutational steps sufficiently to allow mutants to enter the coexitence region, \(P_2\). The derivation of the canonical equation integrates only over those mutants that fall into the invasion-substitution region, \(P_1\). However, since the derivation is done in the limit of small mutational step sizes, this encompasses essentially all of the weight of the mutational kernel, and to good approximation can thus be replaced with the entire \(P_1 \cup P_2\) region. Recall that in the derivation of the canonical equation, deleterious mutations are given zero probability of fixing (from the branching process approximation), and hence this corresponds to integrating over half of the real line \(\mathbb{R}\).
\[\begin{aligned} \frac{\mathrm{d}\bar x}{\mathrm{d}t} = \mu N^*(x) \phantom \cdot \partial_y s(y,x) \Big|_x \int_{P_1} [y-x]^2 M(y-x) \mathrm{d}[y-x] \label{canonical}\end{aligned}\]
The probability of such a mutant occurring in our monomorphic population of trait \(x\) is given by the probability that any mutant occurs (given by the population size, approximately \(K(x)\) at equilibrium, times the birth rate \(b(x)\) times the mutation rate \(\mu\)) times the probability that the mutation lies on the opposite side of the coexistence boundary, which we denote as \(\phi(x)\) and \(\psi(x)\), which is \(\phi(x)\) reflected across the line \(y=x\). This depends on the mutational kernel, \(M(y,x)\) which gives the probability that a mutant arising from a resident population of trait \(x\) bears trait \(y\) that falls withing the coexistance region as follows:
\[\int_{-\infty}^{\phi(x)} M(y,x) \mathrm{d}y + \int_{\psi(x)}^{\infty} M(y,x) \mathrm{d}y\]
Without loss of generality let us assume that the singular strategy \(x^* = 0\) and that we start from some trait value \(x_0 > x^*\). Additionally, let us we assume \(M(y,x)\) is Gaussian in \(y-x\) with mean 0 and variance \(\sigma_{\mu}^2\). However, we do not wish to consider all such mutants \(y\), but only those that survive. We multiply the probability of a mutant having trait \(y\), \(M(y,x)\), times the probability that is survives drift. This probability we quote from the Galton-Watson branching process, \(1-\tfrac{d(y,x)}{b(y,x)}\). Using our model for birth and death rates , the probability of surviving branching is
\[1-\frac{C(y,x)K(x)}{K(y)} := S(y,x) \label{S}\]
which we denote as \(S(y,x)\) as indicated. Hence we modify our integral to consider all mutants \(y < x^*\) which occur with probability \(M(y,x)\) and then survive with probability \(S(y,x)\), and write down the rate \(P(x)\) at which a mutant which leads to branching occurs from our monomorphic population with trait \(x\):
\[P(x) := r \mu K(x)\left( \int_{-\infty}^{\phi(x)} M(y,x) S(y,x) \mathrm{d}y + \int_{\psi(x)}^{\infty} M(y,x) S(y,x) \mathrm{d}y \right) \label{MSerf}\]
Though a little cumbersome, this can be written in nearly closed form (using the error function). For our particular model \(s(y,x)\), we can write the boundaries \(\phi(x)\) and \(\psi(x)\) as
\[\begin{aligned} \phi(x) = x\frac{\frac{\sigma_k^2}{\sigma_c^2}+1}{\frac{\sigma_k^2}{\sigma_c^2}-1} \nonumber \\ \psi(x) = x\frac{\frac{\sigma_k^2}{\sigma_c^2}-1}{\frac{\sigma_k^2}{\sigma_c^2}+1} \label{phipsi}\end{aligned}\]
Then, given a starting condition \(x_0\), we can determine the expected trait \(x\) at time \(t\) of a monomorphic resident population by integrating the canonical equation,
\[\frac{\mathrm{d}x}{\mathrm{d}t} = \frac{1}{2} \mu \sigma_{mu}^2 N^*(x) \partial_y s(y,x)\]
which for our model is
\[\frac{\mathrm{d}x}{\mathrm{d}t} = \frac{-x}{2\sigma_k^2} r \mu \sigma_{\mu}^2 K(x)\]
While no closed form solution exists for \(x(t)\) in this case, if we assume \(x_0 \ll \sigma_k\) then to good approximation we can view \(K(x)\) to be constant over the interval and take our path to be:
\[x(t) = x_0 \exp\left( -r \mu \sigma_{\mu}^2 K_0 t/\sigma_k^2\right)\]
Using this mean path, we can write down an approximation for the rate of a branching mutant occurring as a function of time, \(P(x(t))\). Using this time dependent rate, the probability of not branching after time \(T\) is simply \(\exp\left( - \int_0^T P(t)\mathrm{d}t \right)\). One minus this is the probability branching by time \(T\); i.e. the cumulative density function, hence the probability density function for the waiting time is:
\[\Pi(T) = P(x(T)) \exp\left( -\int_0^T P(x(t)) \right) \label{pdf}\]
and the expected time to complete phases 1 and 2 is \(\int_0^{\infty} T \Pi(T) \mathrm{d}T\). This is our first analytic approximation.