Competition, Trait Variance Dynamics, and the Evolution of a Species’ Range ^††thanks: This work was supported by the NSF grant DMS 1615126 to JRM.

Empirical and theoretical studies have confirmed that competing species often segregate in different parts of the habitat available to them and resist encroachment by competitors in areas where they are established (Pigot and Tobias 2013). By contrast, less evidence for failure to adapt to a continuously varying environment as a range-limiting factor exists (Micheletti and Storfer 2020; Angert et al. 2020; Colautti and Lau 2015; Benning et al. 2019). This is largely due to the difficulty of establishing values for key parameters in evolutionary models, especially the optimum value of a trait under selection as a function of time and space. However, a few empirical studies in both altitudinal and latitudinal environmental gradients support the hypothesis that a population spreading in a heterogeneous environment can be halted by maladaptation, even in the absence of environmental discontinuities (Dawson et al. 2010; Sanford et al. 2006).

Theoretical studies of the joint action of competition and adaptation are therefore valuable, as they can point to conditions under which adaptation, or failure to adapt, become important factors in limiting ranges. Such studies should therefore help researchers develop tests for selection as a range-limiting factor. Ultimately, they should inform environmental management decisions, such as allocating conservation effort among marginal or central populations, planning for the impact of climate change on the structure and survival of natural communities, and controlling invasive species.

The model we present here builds on a sequence of models inspired by an influential idea—genetic swamping—suggested by Haldane (1956) and Mayr (1963). They hypothesized that the intraspecific feedback between core-to-edge maladaptive gene flow and population density of a species can limit its range. Building on the work by Pease et al. (1989), this process was mathematically modeled by Kirkpatrick and Barton (1997) for a single species in a one-dimensional geographic space—as a system of two partial differential equations that represent joint evolution of the species’ population density and the mean value of a quantitative phenotypic trait along a linear environmental gradient in the trait optimum. The quantitative trait is assumed to be directly associated with the fitness of the species’ individuals. Wing loading and body size in flying insects (Takahashi et al. 2016), or specific leaf area and height in plants (Ackerly and Cornwell 2007), are examples of such fitness-related traits. Stabilizing selection on the quantitative trait in the Kirkpatrick and Barton (KB) model enables the species to adapt to new areas and expand its range. However, Kirkpatrick and Barton showed numerically that with sufficiently steep environmental gradients, gene flow from the densely populated center of the species’ range can prevent local adaptation at the periphery and result in an evolutionarily stable range limit—that is, their model demonstrated range pinning through genetic swamping. A derivation of the KB model as well as analytical results on the existence of range expansion traveling waves and some localized stationary solutions for this model are provided by Miller and Zeng (2014), Miller (2019), and Mirrahimi and Raoul (2013).

A number of models have been constructed and (usually numerically) analyzed that build on the KB model. Each of these modeling studies takes some factor or factors omitted from the KB model, checks whether it makes range pinning easier or harder to achieve, and in some cases also studies its effect on quantitative properties like invasion speed. Such variants include individual-based models, models with stochastic differential equations, and models with discrete time and/or space (Duputié et al. 2012; Polechová et al. 2009; Alleaume-Benharira et al. 2006; Polechová 2018; Filin et al. 2008; Bridle and Vines 2007; Polechová and Barton 2015; Kanarek and Webb 2010). Some models incorporate interspecific competition (Leimar et al. 2008; Goldberg and Lande 2006) and some incorporate mutation (Behrman and Kirkpatrick 2011).

Among these varied studies, one key extension of the KB model was developed by Case and Taper (2000). They gave a community context to the KB model by considering multiple species and competitive interactions between them. Through numerical simulations of their model, Case and Taper showed that the interaction of an environmental gradient and gene flow with the frequency-dependent selection generated by phenotypic competition between species can result in range limits at much shallower environmental gradients than without competition. They also showed that their model generates a number of other important evolutionary phenomena, such as species character displacement in sympatry and range shifts in response to climate change. In addition, Case and Taper challenged the conclusions of Kirkpatrick and Barton (1997) by arguing that competition would be an important force setting range boundaries, even in conditions favorable to genetic swamping, for most realistic parameter combinations. Indeed, they argued that realistic trait optimum gradients would rarely be steep enough to limit ranges in the absence of competition (Case and Taper 2000).

Shortly after Case and Taper’s study appeared, Barton extended the original KB model, in which additive genetic variance was held constant, to a model in which this variance was allowed to evolve along with population density and trait mean (Barton 2001). He developed different variants of the model, in which he also incorporated mutation as a constant forcing term in the equations governing the variance. He concluded from numerical simulations that range pinning via genetic swamping was impossible for these models. Instead, he found that genetic variance inevitably rose past a critical value, providing “fuel” that allowed even a temporarily pinned population to spread over the whole domain. In the same book chapter, Barton raised two arguments against Case and Taper’s conclusion that competition would usually play a greater role than genetic swamping in setting range limits. Specifically, he argued that adaptation typically involves several traits, which would increase the selective pressures in a swamping model, and that Case and Taper were too restrictive in their judgment as to how steep a spatial trait optimum gradient might realistically be (Barton 2001).

Juxtaposing the work of Case and Taper (2000) with that of Barton (2001) presents a conflict that must be resolved. That range limits can arise in continuously varying environments is well established (Sexton et al. 2009). When Barton made the KB model (arguably) more realistic by taking account of mutational and other changes in trait variance, the revised models predicted that such range limits could not be set by genetic swamping (Barton 2001). This would imply that other key drivers of range pinning must exist. Competition suggests itself as such a driver, particularly in view of Case and Taper’s findings (Case and Taper 2000). However, Case and Taper held variance constant in their model, so the challenge to the genetic swamping hypothesis posed by Barton’s findings (Barton 2001) cannot be settled merely by appealing to Case and Taper. Rather, we must test whether competition, combined with selection favoring an optimum trait value that varies in space, can still set range limits when trait variance is allowed to evolve together with population density and trait mean.

Here, we resolve the tension between the disparate conclusions of Kirkpatrick and Barton (1997), Case and Taper (2000), and Barton (2001) by incorporating both competition and nonconstant trait variance in a model of adaptation and spread in an environmental gradient. We take particular care to establish plausible ranges for parameter values, since these ranges play a role in the debate over whether genetic swamping is commonly a dominant influence on range limits. We explore the behavior of solutions to our model numerically, finding that the results of Case and Taper (2000) are robust when the assumption of constant trait variance is eliminated and even when mutation, as a perpetual source of genetic variance, is incorporated in the model through a forcing term as in the work of Barton (2001). We focus primarily on the establishment of stable boundaries between the ranges of two competitors initially established in allopatry. However, we also briefly investigate other scenarios, including the ability of a newly arrived species to establish itself and spread in a habitat already occupied by a competitor.

Biological species are diverse in their level of abundance, growth rate, dispersal range, and the nature of their functional traits. This implies that an appropriate choice of units for the physical quantities of the model such as time, space, populations density, and trait values must be species-dependent. For this, we first select one of the species from the community of $\mathrm{N}$ species in the model as a representative species. The physical units are then chosen (in principle) based on specific measurements made on this species. The representative species can be selected, for example, as the species that is best adapted to the environment, is most widely spread, or has the widest trait niche among the community.

Let $\mathtt{T}$ denote the unit of time. We set $1\;\mathtt{T}$ to be equal to the mean generation time of the representative species. This is a natural choice of time unit to analyze population dynamics of species over evolutionary time scales; see, for example, the time unit used by Estes and Arnold (2007). It also makes the model compatible with common experimental approaches in estimating parameters such as the the strength of phenotypic selection, which is often estimated by measuring changes in the mean phenotype of the population in one generation. Moreover, choosing generation time as unit of time can make the predictions of the model comparable with the results obtained from discrete-time individual-based models which usually describe the evolution of the population at generation time steps; see, for example, the models used by Polechová and Barton (2015) and Bridle et al. (2010).

To set the unit of space, denoted by $\mathtt{X}$, we first consider a one-dimensional habitat, that is, $x\in\mathbb{R}$. It can be seen in equations (1)-(3) that rescaling the space as $x\mapsto kx$ leaves the equations unchanged, provided the diffusion coefficients $\mathrm{D}_{i}$ are rescaled accordingly as $\mathrm{D}_{i}\mapsto\mathrm{D}_{i}/k^{2}$. Using this flexibility in the equations and having set the unit of time, we choose the unit of space so that the dispersal coefficient of the representative population becomes unity, that is, $1\;\mathtt{X}$ is the root mean square dispersal distance of the population in $1\;\mathtt{T}$ divided by $\sqrt{2}$. Since the population dispersal may vary at different locations of the habitat, the measurements for setting the unit can be done based on a local subpopulation at the core of the population or at regions where dispersal is not affected by environmental barriers. For multi-dimensional habitats, the same approach can be used to set the unit of space for each spatial dimension independently. We denote the units associated with all dimensions by $\mathtt{X}$, noting that $\mathtt{X}$ may refer to different physical scales for different spatial axes.

Similarly, rescaling the population density of the species as $n_{i}\mapsto kn_{i}$ does not change the equations of the model provided the carrying capacity of the environment is rescaled accordingly as $\mathrm{K}_{i}\mapsto k\mathrm{K}_{i}$. Therefore, having set the unit of space, we choose the unit of measurement for population abundances, denoted by $\mathtt{N}$, so that the carrying capacity of the environment for the representative population becomes unity. That is, $1\;\mathtt{N}$ is equal to the carrying capacity of the environment for $1\;\mathtt{X}^{\mathrm{m}}$ unit of habitat volume. The required measurement of the carrying capacity can be done locally at the core of the representative population where it has the largest and most stable population density.

Finally, we denote the unit of measurement for the quantitative trait by $\mathtt{Q}$, and set $1\;\mathtt{Q}$ to be equal to one standard deviation of the trait values at the core of the representative population, where the population likely shows highest variance in the individual’s trait values. This is a common choice of unit for quantitative traits, which provides generality for quantitative models of evolutionary processes by making them independent of the diverse nature of the quantitative traits across different species (Kingsolver et al. 2001; Lande and Arnold 1983; Estes and Arnold 2007; Kirkpatrick and Barton 1997; Case and Taper 2000).

Based on the units chosen in Section 3.1 for measuring physical quantities of the model, plausible ranges of parameter values can be suggested as follows.

The amount of increase in genetic variance of phenotypic traits per one generation of mutation—after being standardized with (divided by) the estimates of environmental variance of the trait—is known as mutational heritability in the literature. Estimates of this dimensionless quantity are available for a variety of species (Houle et al. 1996) and can be used here to obtain biologically reasonable values for $\mathrm{U}$. Note that the environmental trait variance is $(1-\mathrm{H}^{2})$ times the total phenotypic variance, where $\mathrm{H}^{2}$ denotes the broad sense heritability (Visscher et al. 2008). Moreover, since the standard deviation of the phenotypes at the core of a representative population is chosen here as the unit of trait, as described in Section 3.1, the typical value of the total trait variance in our model is expected to be approximately $1\,\mathtt{Q}^{2}$. Therefore, the estimates of mutational heritabilities provided by Houle et al. (1996, Table 1), after being multiplied by the factor $(1-\mathrm{H}^{2})$, give plausible values for $\mathrm{U}$.

The standardized values given by Houle et al. (1996, Table 1) are approximately in the range $[0.1,30]\times 10^{-3}$. Typical values of heritability for fitness-related traits, as given by Visscher et al. (2008, Figure 1) for a number of different species, range from $0.05$ to $0.3$. These results can suggest $[0.1\times(1-0.3),30\times(1-0.05)]\times 10^{-3}\approx[7\times 10^{-5},0.029]\,\mathtt{Q}^{2}/\mathtt{T}$ as an approximate range of values for $\mathrm{U}$. Specifically, the standardized estimates of mutational rate of increase provided by Houle et al. (1996, Table 1) for Daphnia pulex are in the range $[0.8,3.3]\times 10^{-3}$, and for Drosophila melanogaster are in the range $[0.1,13.5]\times 10^{-3}$. Values of heritability given by Visscher et al. (2008, Figure 1) for fitness-related traits in Daphnia and Drosophila are approximately $0.3$ and $0.18$, respectively. Noting that not all the traits considered by Houle et al. (1996, Table 1) are fitness traits, these results give the estimates $[0.8\times(1-0.3),3.3\times(1-0.3)]\times 10^{-3}\approx[0.0006,0.0023]\,\mathtt{Q}^{2}/\mathtt{T}$ and $[0.1\times(1-0.18),13.5\times(1-0.18)]\times 10^{-3}\approx[8\times 10^{-5},0.01]\,\mathtt{Q}^{2}/\mathtt{T}$ for the range of values that $\mathrm{U}$ can take for Daphnia and Drosophila, respectively.

In addition to providing estimates of the mutational rate of increase in trait variance relative to the environmental variance, Houle et al. (1996) also provide estimates of the rate relative to the standing genetic variance. Noting that genetic variance is equal to $\mathrm{H}^{2}$ times the total phenotypic variance (Visscher et al. 2008), we can use similar calculations as given above to obtain additional estimates of values for $\mathrm{U}$. The estimates provided by Houle et al. (1996, Figure 4) are in the range $[10^{-2.75},10^{-1.25}]\approx[0.002,0.056]$, relative to the genetic variance. Therefore, considering the values of heritability from $0.05$ to $0.3$, as given above, we obtain a range of possible values for $\mathrm{U}$ as $[0.002\times 0.05,0.056\times 0.3]\approx[0.0001,0.017]\,\mathtt{Q}^{2}/\mathtt{T}$. It should be noted that, since in our model we implicitly assume that variation in phenotypes is mainly caused by genetic effects, that is $\mathrm{H}^{2}\approx 1$, we expect to use slightly larger values than those estimated here for te parameter $\mathrm{U}$ in our model.

An extremely rough estimate for the range of values of $\mathrm{U}$ can be obtained from the spatially homogeneous equilibrium value of trait variance in (1)–(3). For a perfectly adapted solitary species at spatially homogeneous equilibrium, and in the absence of environmental gradients and intraspecific competition, the equilibrium trait variance $v^{\ast}$ is maintained by the mutation-selection balance $\mathrm{U}=\mathrm{S}{v^{\ast}}^{2}$; see equation (14) and the discussions in Section 4. With a typical value of $v^{\ast}=1\,\mathtt{Q}^{2}$, this gives estimates of values for $\mathrm{U}$ equal to the strength of stabilizing selection $\mathrm{S}$, which as discussed above may range from $0$ to $2$. However, gene flow over an environmental cline and intraspecific competition are indeed prominent sources of producing genetic variation, which were ignored in obtaining this rough estimate of the range of values for $\mathrm{U}$. Therefore, a value of $\mathrm{U}=2\,\mathtt{Q}^{2}/\mathtt{T}$ will be too extreme to be realistic, although in shallow environmental gradients we may expect the value of $\mathrm{U}$ to be relatively close to $\mathrm{S}$.

Finally, considering altogether the estimates and sample values provided above, we suggest a plausible range of values for $\mathrm{U}$ as $[0,0.2]\,\mathtt{Q}^{2}/\mathtt{T}$, with a typical value of $0.02$. However, we set $\mathrm{U}=0$ in all of the computational studies presented in this paper, except for the results discussed in Remark 4.3 in Section 4. This is because the inflation caused by mutation in species’ trait variance does not affect the conclusions of our numerical studies; see Remark 4.3 for details.

We assume the environmental trait optimum has a constant gradient $\mathrm{d}_{x}\mathrm{Q}=0.2\;\mathtt{Q}/\mathtt{X}$, which is the slope of the black line in Figure 2(b). The species is introduced at the center of the habitat with a density given as $n(x,0)=0.5\operatorname{sech}(|x|/\sqrt{2})$. It is assumed that this initial population is perfectly adapted to the environment at the center, that is, $q(0,0)=\mathrm{Q}(0)$, and has a linearly varying trait mean of slope $\partial_{x}q(x,0)=0.6\,\mathrm{d}_{x}\mathrm{Q}$ for all $x\in\Omega$. We further assume that the initial population has a constant trait variance of $v(x,0)=1\;\mathtt{Q}^{2}$.

Figure 2 shows the solutions of the model (12)–(14) over the computation time horizon of $T=50\;\mathtt{T}$. It can be seen that the species’ population density initially grows to an upper limit relatively fast, and then the population invades the entire habitat in the form of a traveling wave. As the population spreads over the habitat, it successfully adapts to new areas in response to the force of natural selection, and its mean trait gradually converges to the optimal trait. The initially constant profile of the population’s trait variance evolves quickly to a bell-shaped profile, showing a larger trait variance at the core of the population with a gradual decline towards the edges. The maximum trait variance at the population center evolves fairly slowly to an upper bound as the population expands its range across the habitat.

To investigate further the invasive range dynamics of this solitary species, a sample curve is highlighted in the solutions shown in Figure 2. At the effective edges of the population, that is considered to correspond to the inflection points on the curve of the population density, the diffusion term $\mathrm{D}\partial_{x}^{2}n$ in (12) is zero. However, the mean intrinsic growth rate $G$, given by the term inside the parentheses in (12) and shown in Figure 2(d), is positive at these edge positions. This implies that the population density $n$ is not stationary at the edges, which confirms the fact that the sample solution curve is a traveling wave.

The curve of trait mean highlighted in Figure 2(b) shows substantial adaptation at the core of the population, but considerable failure in adaptation near the edges. This adaptation profile is mainly due to the effect of gene flow and is observed, after few generations, even if the initial population is perfectly adapted everywhere, that is, even if $q(x,0)=\mathrm{Q}(x)$ for all $x\in\Omega$. This profile can be explained as follows. The density of the population is nearly uniform at its core, and hence there is a fairly symmetric gene flow to a core location from its adjacent areas located on the upper and lower parts of the cline. Due to this symmetry, gene flow does not significantly affect the mean value of the trait at core locations and adaptation is successfully maintained by natural selection. Near the edge, however, the population density varies sharply and the gene flow is highly asymmetric. At the left edge shown in Figure 2(b), for example, gene flow is predominantly from the upper points on the cline. Therefore, near the left edge, the mean value of the trait is increased above the optimal value due to the effect of this asymmetric gene flow. As a result, the curve of trait mean is gradually flattened and departs from the optimal curve as it approaches the left edge. A similar process, but in opposite direction, flattens the curve below the optimal value at the right edge.

The bell shape of the trait variance is a result of gene flow and the specific adaptation pattern described above. At core locations, the population is well-adapted to the optimal cline and, due to the relatively large gradient of the cline, gene flow from adjacent areas generate large phenotypic variations among central individuals. Near the edges, however, the population fails to adapt to the optimal cline and the curve of trait mean flattens. Therefore, due to the low gradient in the trait mean near the edges, gene flow from adjacent areas does not substantially contribute to phenotypic variations in marginal individuals. As a result, trait variance constantly declines from core to edges, in parallel with the decline in the gradient of the trait mean.

The profile of variations in the trait mean and the trait variance described above are consistent with experimental observations available in the literature. For example, Takahashi et al. (2016) have shown geographic variations in abdomen length and wing loading of two closely related species of damselflies along a latitudinal cline in Japan. The patterns of variation in both of these phenotypic traits indicate that the species are well adapted to the environmental gradient at the core of their range, whereas they are significantly maladapted at their range margins. Additionally, variations in heterozygosity and the Garza-Williamson index, as the two indicators of genetic diversity measured by Takahashi et al. (2016), show drastic decline in species’ genetic and phenotypic variation at their range margins where maladaptation is observed.

A major difference between the predictions of the present model and predictions of the models that assume constant trait variance (Kirkpatrick and Barton 1997; Case and Taper 2000) appears in response to steep environmental gradients. To show this, we first repeat the computation of Section 4.1 for a species under a steep environmental gradient of $\mathrm{d}_{x}\mathrm{Q}=2\;\mathtt{Q}/\mathtt{X}$, that is $10$ times larger than the typical gradient considered in Section 4.1. Here, we initialize the computation with a better adapted population with $\partial_{x}q(x,0)=0.9\,\mathrm{d}_{x}\mathrm{Q}$, to avoid the numerical singularities that would otherwise occur at initial iterations of the computation—mainly due to large trait mismatch, $|q-\mathrm{Q}|$, occurring near the boundary of the habitat because the environmental gradient is very steep. The rest of the computation parameters take the same values as used in Section 4.1.

The population density and trait variance of the species is shown in Figure 3. It can be seen that the species is still able to expand its range in the form of a traveling wave, but with a lower speed compared to the wave speed at the typical gradient value. More importantly, the initially small value of the trait variance, $v(x,0)=1\,\mathtt{Q}^{2}$, evolves quickly to a significantly larger value. In fact, since the population here is adapting itself to an optimal phenotype with much steeper variation along the habitat, the mean phenotypes diffused to a given location within the species’ range have a much wider range of values. As a result, gene flow in this case generates large phenotypic variations in the population. Similar to the variance profile described in Section 4.1, the trait variance declines at the range margins due to the effect of asymmetric gene flow. Note that, for $\mathrm{d}_{x}\mathrm{Q}=2\,\mathtt{Q}/\mathtt{X}$ considered here, the models that assume constant trait variance show an evolutionarily stable limited range, for which the mean intrinsic growth rate vanishes to zero at the inflection points on the range edges, and remains negative beyond. Examples of such limited range dynamics are provided by Kirkpatrick and Barton (1997) and hence are not presented here.

Next, to see in more detail the differences resulting from the evolution of phenotypic variance in the range dynamics of the species, we additionally generate a constant-variance version of the model and repeat the computations described above for both models, and with different values of the environmental gradient. The constant-variance model is obtained by fixing $v(x,t)=\mathrm{V}_{\textsc{p}}$ in (12)–(14), that is, by omitting (14) and replacing $v$ by $\mathrm{V}_{\textsc{p}}$ in (12) and (13). The resulting model is equivalent to the model presented by Case and Taper (2000) and Kirkpatrick and Barton (1997). We set $v(x,0)=1\,\mathtt{Q}^{2}$ for our variable-variance model and, correspondingly, $\mathrm{V}_{\textsc{p}}=1\,\mathtt{Q}^{2}$ for the constant-variance model.

For each value of the environmental gradient, solutions of the equations are calculated for a sufficiently large period of time. The approximate speed and amplitude³³3Here, we define the amplitude of a population’s traveling wave solution as the (steady) peak value of the population’s density during its range expansion regime. In the results presented in Figure 4(b), the traveling wave amplitudes are approximately calculated as the density of the population at its center when it has reached a nearly steady value. of the population’s traveling waves are then shown in Figure 4. Note that, at every gradient value, the position of a point on the edges of the population density waves initially changes nonlinearly in time. However, when the effect of initial transient dynamics vanishes, the variations in the edge positions eventually become linear with respect to time. The slope of these linear variations is shown in Figure 4(a) as an approximate value of the wave speed at each gradient value. Moreover, the amplitude of the traveling wave at the final time of computations for each gradient value is shown in Figure 4(b) as an approximate value of the wave amplitude.

Figure 4 shows that the constant-variance model predicts a relatively sharp decline in the speed of invasion waves as the environmental gradient increases, but no change in waves’ amplitude. Kirkpatrick and Barton (1997) and Case and Taper (2000) show that a stable limited range is formed if the gradient is increased beyond the value at which wave speed becomes zero. Eventually, the population becomes extinct if the gradient is further increased to very large values. However, removing the assumption of constant phenotypic variance, as in the present model, results in essentially different predictions. It can be seen in Figure 4 that, although the speed of invasion waves declines (almost linearly) when the environmental gradient increases, the species is still able to expand its range under much steeper gradients than what is predicted under the constant variance assumption. Nonetheless, the amplitude of the species’ expansion waves declines substantially at steep gradients and the population eventually vanishes at extreme gradients. This implies that physical barriers, such as mountains or oceans, which impose extreme environmental gradients on the species’ habitat can still prevent the species’ range expansion.

The approximate wave amplitudes shown in Figure 4(b) can be alternatively calculated using the equations of a fully adapted homogeneous population. As the range dynamics of the species shown in Figures 2 and 3 suggests, population density and trait variance eventually become homogeneous over $\Omega$ as $t\to\infty$, and trait mean converges uniformly to $\mathrm{Q}$. We denote this equilibrium state by $(n^{\ast},q^{\ast},v^{\ast})$, where $n^{\ast}>0$ and $v^{\ast}\geq 0$. At this equilibrium, we have $\partial_{x}n^{\ast}=0$, $q^{\ast}=\mathrm{Q}$, and $\partial_{x}v^{\ast}=0$. Moreover, $\partial_{x}^{2}q^{\ast}=0$ since $\mathrm{Q}$ is assumed to be linear here. Therefore, from (12)–(14) we obtain

and

Letting $v^{\ast}=\mathrm{V}_{\textsc{p}}=1$ in (15) for the constant-variance model, we simply obtain $n^{\ast}=1.062$ as the constant wave amplitude shown in Figure 4(b). For the variable variance model, the cubic algebraic equation (16) has one and only one positive root for all nonzero values of $\mathrm{U}+2\mathrm{D}|\mathrm{d}_{x}\mathrm{Q}|^{2}$. The graph of this nonzero root with respect to changes in the gradient $\mathrm{d}_{x}\mathrm{Q}$ appears to be approximately a straight line with positive slope. Moreover, substituting this root into (15) for different values of $\mathrm{d}_{x}\mathrm{Q}$ provides a graph of $n^{\ast}$ with respect to $\mathrm{d}_{x}\mathrm{Q}$. With the parameter values considered in this section, this graph gives an approximate curve of wave amplitudes as shown in Figure 4(b). Note that $n^{\ast}=0$ in (15) when $v^{\ast}=2\mathrm{R}/\mathrm{S}$, which is a solution of (16) with $|\mathrm{d}_{x}\mathrm{Q}|=\sqrt{(2\mathrm{R}^{2}/\mathrm{S}\mathrm{D})-\mathrm{U}/2\mathrm{D}}=6.325$. This gives the value of the environmental gradient at which the wave amplitude becomes zero and the population fails to survive.

An abrupt change in the curve of optimal phenotype may occur due to a rapid climate change and can largely affect a species’ geographic range and abundance, particularly when it occurs frequently. To observe the predictions of the model in response to these changes in climate, here we initialize the computations at $t=0\,\mathtt{T}$ using the solution curves computed in Section 4.1 at $t=4\,\mathtt{T}$. However, we uniformly shift up the line of trait optimum $\mathrm{Q}$ by a factor of $5\,\mathtt{Q}$ at $t=0\,\mathtt{T}$. As a result, a high level of maladaptation is immediately induced in the population and the population’s density and expansion speed decline quickly. These impacts are more severe near the right margin of the population’s range, where, as discussed in Section 4.1, the population initially has a trait mean below the trait optimum and hence faces a greater maladaptation as the optimum is shifted up.

The impact of the climate change described above is transient and gradually dissolves as the species adapts itself to the new environmental gradient. However, the impacts of climate change can last longer and become more severe if rapid changes in climate occur more frequently. For instance, Figure 5 shows the results obtained when abrupt fluctuations of magnitude $5\,\mathtt{Q}$ occur periodically in the trait optimum. If the period of these fluctuations is sufficiently small, as in Figure 5, the abundance and expansion speed of the species is steadily affected by the fluctuations. The species’ range margins, particularly the right margin, advance slower under the climate fluctuations and the species’ density remains significantly below its environmental carrying capacity. Moreover, note that the pattern of fluctuations will shape the longterm profile of the population’s trait mean. Here, fluctuations are in the form of a square wave, with $\mathrm{Q}$ shifting up by a factor of $5\,\mathtt{Q}$ for the first half of the fluctuations’ period, and then shifting back to its initial value for the second half. Regulated by the evolutionary force of natural selection, the trait mean then converges to a line which is above the initial $\mathrm{Q}$ by a factor equal to half of the periodic shift, that is, $5/2\,\mathtt{Q}$. This can be seen in Figure 5(b).

The impacts of the periodic fluctuations described above can be more severe if the fluctuations have larger amplitudes or the species evolves under stronger phenotypic selections. For instance, if we repeat the computations of Figure 5 with a stronger selection, $\mathrm{S}=0.6\,{\mathtt{Q}}^{-2}/{\mathtt{T}}$, we obtain the results shown in Figure 6. It can be seen that, in this case, the population fails to withstand the climate fluctuations and becomes extinct. Figure 6(b) shows that the natural selection is still regulating the trait mean at the best possible value, that is $5/2\,\mathtt{Q}$ above the initial $\mathrm{Q}$. But, this does note give the population the fitness required for survival under such a strong selection.

It is shown by Case and Taper (2000) that interspecific competition can effectively limit the range of two competing species that would expand indefinitely in the absence of competition. Here, we show that this result is not significantly affected by the evolution of the intraspecific phenotypic variance, and similar predictions still hold when the constant variance assumption is removed. For this, we consider two populations of species that are initially distributed allopatrically and are both perfectly adapted to the environment at their center. Similar to the initial values considered for the single species of Section 4.1, we set $\partial_{x}q_{1}(x,0)=\partial_{x}q_{2}(x,0)=0.6\,\mathrm{d}_{x}\mathrm{Q}$ and $v_{1}(x,0)=v_{2}(x,0)=1\;\mathtt{Q}^{2}$ for all $x\in\Omega$.

Figure 7 shows the solutions of the model for $T=50\,\mathtt{T}$. Each species expands its range to the edge of the habitat on the side where the species initially resides. However, in the middle of the habitat where the two species meet, they prevent each other’s progress. At the beginning, when the species are apart from each other, their range dynamics is basically the same as the dynamics of a solitary species shown in Section 4.1, that is, their trait mean gradually converges to the trait optimum as the species adapt and advance to new areas. Therefore, when the species first meet at the middle of the habitat, the individuals of either species that are near the interface of the two populations have relatively close phenotypes. As a result, a strong interspecific competition is initiated between these individuals, according to the competition kernel specified by equation (23) in Appendix A.4. This competition decreases the fitness of the populations at their interface. Hence, the density of the populations declines over their interface, and this in turn intensifies the effect of asymmetric gene flow within each population at the interface. By the same process as described for a single species in Section 4.1, this asymmetric gene flow gradually flattens the trait mean curve of each population over the interface, so that the curves depart from the trait optimum curve in opposite directions. This can be seen through the highlighted curves in Figure 7(b).

The interaction described above between gene flow and interspecific competition continues until the level of maladapted phenotypes in each species’ peripheral populations inside the interface becomes so extreme that it prevents local adaptation and hence stops species’ range expansion through the interface. Consequently, a region of sympatry is formed between the species at the middle of the habitat, over which the species’ population density monotonically declines to zero. Figure 7(a) shows the formation of the associated range limits. Moreover, as shown in Figure 7(d), the species exhibit significant character displacement in sympatry. These observations are consistent in general with those given by Case and Taper (2000). Note that smaller values of the variance of individuals’ phenotypic utilization, $\mathrm{V}_{i}$, restrict the interspecific competition to individuals with closer phenotypes. As a result, the overall interspecific competition becomes weaker and the region of sympatry becomes wider. The width of the region of sympatry increases also with shallower environmental gradients, which result in less extreme asymmetry in the gene flow. If the environmental gradient is zero, the two species of Figure 7 eventually become sympatric over the entire available habitat.

The highlighted curve in Figure 7(c) shows that the phenotypic variance within each species evolves to a bell-shaped curve. This is a consequence of asymmetric gene flow at both edges of the species’ range, as explained for the single species of Section 4.1. Note that the large peaks in the curves shown in Figure 7(c) occur where the population has an infinitesimal density. Such regions are not practically considered within the range of the species, and values of trait mean and variance over these regions are of no biological meaning.

Generically, the region of sympatry originated by interspecific competition, as described above in Section 5.1, does not remain stationarily centered between the two species. The stationary region that is seen in Figure 7(a) is simply due to the identical choice of parameter values for both species. Any differences in the parameters that change the competition balance between the two species cause the interface to constantly move towards the competitively weaker species. If the imbalance is large enough, the weaker species goes extinct after being ultimately pushed to the boundary of the habitat. However, the dynamics of this competitive exclusion process is relatively slow. For instance, if we repeat the computations of Figure 7, but this time with $\mathrm{R}_{1}=3\,\mathtt{T}^{-1}$ and $\mathrm{R}_{2}=2\,\mathtt{T}^{-1}$, we see that a region of sympatry is formed quickly in the middle of the habitat within $10\,\mathtt{T}$, but it then moves rather slowly towards the right boundary. The second species reaches the right boundary and begins to decline in density approximately at $t=270\,\mathtt{T}$. It practically becomes extinct at about $t=460\,\mathtt{T}$. This evolutionary dynamics can indeed be much slower if the difference between $\mathrm{R}_{1}$ and $\mathrm{R}_{2}$ is made smaller.

A similar, but biologically more interesting competitive exclusion occurs when an invasive species is introduced to a habitat that is already occupied by a well-adapted species. To see this, we assume that the final population of the solitary species of Section 4.1 at $t=50\,\mathtt{T}$ is occupying the habitat here at $t=0\,\mathtt{T}$, labeled as the $1$st species. Moreover, we assume a faster growing $2$nd species with $\mathrm{R}_{2}=2.5\,\mathtt{T}^{-1}$ is introduced at $t=0\,\mathtt{T}$ at the center of the habitat, with the same initial conditions used for the species of Section 4.1. With these initial populations, we compute the solutions over the time horizon of $T=500\,\mathtt{T}$. The results are shown in Figure 8. It can be seen that, after a transient initial reduction in its density, at about $t=20\,\mathtt{T}$ the introduced species starts growing constantly and expanding its range. The established species, which is significantly weaker than the invasive species due to its smaller maximum growth rate, declines constantly in density over the expanding range of the invasive species. When this species’ central population at the origin becomes extinct, two separate regions of sympatry are formed on opposite sides of the habitat. From this point on, these regions move in opposite directions towards the boundary of the habitat, and the introduced species eventually establishes itself over the entire habitat by fully excluding the preexisting species. This eco-evolutionary process, however, occurs very slowly.

If the two species of Section 5.1 differ from each other, but the difference between them is rather minor, then a stable equilibrium can exist in which the weaker species is not entirely excluded from the habitat, but it appears only over a limited extent adjacent to the boundary. To see this, we consider two populations of species which differ only in their maximum growth rate, with $\mathrm{R}_{1}=2.1\,\mathtt{T}^{-1}$ and $\mathrm{R}_{2}=2\,\mathtt{T}^{-1}$. We initialize a simulation with these two populations being allopatrically distributed and perfectly adapted to the environment everywhere, that is, $\partial_{x}q_{1}(x,0)=\partial_{x}q_{2}(x,0)=\mathrm{d}_{x}\mathrm{Q}$. Moreover, we set $v_{1}(x,0)=v_{2}(x,0)=1\;\mathtt{Q}^{2}$ for all $x\in\Omega$. Figure 9 shows the solutions of the model with these initial populations over a long time horizon of $T=5000\,\mathtt{T}$.