Biological advection and cross-diffusion with parameter regimes

1.   Introduction

A realistic dispersal of biological organisms should take account the interaction between individuals and the response to the environmental variations. The importance of having such a diffusion model has been emphasized by many researchers (see Skellam ^[25,26] and Okubo & Levin [24, Chapter 5]). In a reaction-diffusion equation, such interactions and responses are given in the reaction term and hence it is important to appropriately design the dispersal part depending on the reaction term. If dispersal and reaction parts are not related appropriately, the resulting system may give us a wrong conclusion even if each of them are realistic and correct. The starvation driven diffusion (or SDD for brevity) is such a dispersal model that includes the both components in a simplest way by simply increasing dispersal rate if starvation starts (see ^[3,15,16]). By choosing an appropriate starvation measure for a given reaction term, we may obtain an appropriate relation between the two terms.

The starvation driven diffusion can be called a Fokker-Plank type diffusion or, more precisely, an asymmetric diffusion and is written by

where $\gamma$ is the migration rate or motility (see discussions in Section 2). The purpose of this paper is to show how to derive biological advection and cross-diffusion models from such an asymmetric diffusion operator. We will eventually obtain a systematic method to derive advection or cross-diffusion terms which exploit the population dynamics of a given problem. This derivation enables us to restrict meaningful parameter regimes of advection and cross-diffusion models in the context of biological random dispersal. If parameters are beyond these regimes, the advection or cross-diffusion terms can be justified only by a randomness breaking mechanism.

Let $u(x, t)$ be the population density at position $x\in\Omega$ and time $t > 0$. Let $\gamma$ be the migration rate of the species to depart its habitat. We assume that the migration rate $\gamma$ depends on the amount of resource $m$, the population of other species $v$, and of its own species $u$, i.e, $\gamma = \gamma(u, v, m)$. If the direction of migration is randomly chosen, the corresponding dispersal model is given by an asymmetric diffusion equation,

The four terms in the right are regular diffusion, self-diffusion, cross-diffusion, and advection in the order. We will obtain simplified advection, self-diffusion, and cross-diffusion models from this nonlinear diffusion equation by taking certain coefficients to be constant for simplification. A different choice ends up in a different dispersal with different behavior.

The hypothesis in this derivation is that biological organisms start to migrate if their habitat turns into a harsh place. Various elements can be involved in the migration of real organisms. However, the only thing given in a reaction-diffusion equation is the population dynamics and we can measure the environmental harshness only from the reaction term.^* The specific growth rate is an indicator that shows how good the habitat is. For example, let the population density $u(x, t)$ satisfy a reaction diffusion equation,

^*In this context, we may end up with a wrong dispersal model if the given population dynamics is unreasonable. Hence, the derivation can give a clue about the validity of a reaction term.

where $f$ is the population dynamics of the species. We introduce two ways to set up the migration rate $\gamma$. First, we may assume that $\gamma = \gamma(\tau)$ is a function of the specific growth rate $\tau$ defined by

The reaction term $f$ is the growth rate of the whole colony level and $\tau$ is of the individual level. This specific growth rate is the quantity that matters to each individual. If the specific growth rate $\tau$ increases, each individual has a better chance to grow and hence wants to stay. Therefore, we assume $\gamma$ is a decreasing function of $\tau$, i.e., ${\partial\over \partial\tau}\gamma\le0$. Second, we may assume $\gamma = \gamma(s)$ is an increasing function of a starvation measure $s > 0$. The choice of a starvation measure depends on the population dynamics $f$ and we consider three examples in this paper. We introduce the driven dispersal models in follows as a summary.

The first example is the logistic equation,

where $m = m(x) > 0$ is the carrying capacity at $x\in\Omega$. We are interested in spatial heterogeneity in $m$ and hence $m = m(x)$ is nonconstant function in this paper. The specific growth rate is written as

where $s: = {u\over m}$ is the starvation measure of this example and $\gamma$ is an increasing function of $s$. In Section 3, we derive an advection model,

after a simplification of (1.2). Notice two things. First, this is a logarithmic model which appears naturally throughout the paper. Second, the advection coefficient $\chi_0$, which is often called a chemosensitivity, has a parameter regime $0 < \chi_0 < 1$. Such a parameter regime does not exist in classical chemotaxis theory (see Appendix A).

A variation of the logistic population dynamics is often used as in

In this case, the population dynamics is written as $u(m-u) = mu(1-{u\over m})$, where $m(x)$ is the growth rate and the carrying capacity at the same time. This unusual situation ends up with an interesting dispersal model. We derive two advection models,

Find that the advection term in (1.4) is independent of $u$. Since this advection term is valid only when $u > 0$, the characteristic function ${\bf 1}_{\{u > 0\}}$ is multiplied. The coefficient $\chi_0$ is similarly bounded by one. On the other hand, the other model (1.5) contains usual self-diffusion and advection.

The second example is the Lotka-Volterra equations for two competing species,

where $v$ is the population density of the other species. The coefficients $b_{ij} > 0$ are inter-species competition rates. The specific growth rate of the first equation is $\tau = {f\over u} = 1-{u+b_{12}v\over m} = :1-s_u$, where $s_u: = {u+b_{12}v\over m}$ plays the role of the starvation measure of the first species and $\gamma$ is an increasing function of $s_u$. In Section 4, we derive three cross-diffusion models,

where the parameter regimes are $0 < \chi_0 < 1$ and $p_0 > 0$. Find that the cross-competition coefficient $b_{12}$ appears in cross-diffusion terms. The first equation (1.6) contains a cross-diffusion term without self-diffusion. The coefficient $\chi_0$ appears in both advection and cross-diffusion terms. The second equation (1.7) is with self-diffusion. The coefficient $p_0 > 0$ has no upper bound and appears in three places of advection, self-diffusion, and cross-diffusion. The third equation (1.8) is a degenerate case.

The third example is the Lotka-Volterra predator-prey equations,

where $w$ is the population density of a predator species. In this model, the starvation measure is not appropriate, but the specific growth rate still is. The specific growth rates for the two species are respectively $\tau_u = a-bw$ and $\tau_w = cu-d$, and the departing probabilities $\gamma = \gamma(\tau_u)$ and $\gamma = \gamma(\tau_w)$ are decreasing functions of $\tau_u$ and $\tau_w$. In Section 5, we obtain cross-diffusion terms for predator and prey species,

Find that the coefficients, $a, b, c$, and $d$, do not appear in cross-diffusion terms. The prey has a cross-diffusion term with positive sign in (1.9) and hence the cross-diffusion gives a flux moving away from predator. However, the cross-diffusion coefficient of the predator species in (1.10) is negative and gives a flux moving toward prey. These signs are automatically decided by the signs of specific growth rates. In either case, there is no self-diffusion term because there is to self-competition in the population dynamics.

The migration of biological species is a survival strategy when the environment becomes harsh. There have been intensive studies on the evolution of dispersal strategies of biological organisms. The basic dispersal model is the one with constant diffusivity and only the size effect of constant diffusivity is documented well. Mathematical analysis and numerical simulation show that a smaller dispersal rate is selected if the environment is spatially heterogeneous and that a larger dispersal rate can be selected if there is a temporal fluctuation of the environment (see ^{[1,4,6,7,9,10,11,12,20,22,23,24,27]}). However, biological organisms change their dispersal rate if the environment is changed. The essence of biological dispersal is in such adaptation and the starvation driven diffusion model (1.1) contains it. The advection, self-diffusion, and cross-diffusion models obtained in this paper are from this adaptability.

2.   Random diffusion versus symmetric diffusion

Diffusion is often modeled by the heat equation. Let $d > 0$ be the diffusivity. If $d$ is constant, the heat equation can be written as

which are all equivalent. However, if the diffusivity $d$ is not constant, they are all different. Finding the correct diffusion law in a spatially heterogeneous environment is an unresolved task. Remember that the diffusivity $d$ is decided by many microscopic level components such as mean free path, jumping time, particle speed, turning frequency, etc. Recently, Kim and Seo ^[17] showed that diffusivity alone is not enough to explain heterogeneous diffusion and derived a new diffusion law,

where $\mu$ is the turning frequency.

The focus of the dispersal model in this paper is in another substantial component for the dispersal of biological organisms which is the migration rate. Each organism may choose to leave or stay at their habitat. The diffusion equation (1.1) is about such biological component. In this section, we discuss the meaning of the diffusion operators in this context. We will see that Fick's law is not a random dispersal, but asymmetric dispersal (1.1) is.

Consider a linearly ordered patch system given in Figure 1. Let $u_i(t)$ be the population at patch $i$ and time $t > 0$ and $c_{ij}\ (\equiv c_{i\leftarrow j})$ be the migration rate from patch $j$ to patch $i$. We assume $c_{ij}\ne0$ only for $j = i\pm1$. In other words, the population can migration only to one of two adjacent patches. Then, the population $u_i$ is modeled by

where the first two positive terms are incoming flux from the two neighbor patches and the next two negative ones are outgoing flux.

Symmetric diffusion. The dispersal is called symmetric if $c_{ij} = c_{ji}$. This is a case that the migration rate from patch $i$ to patch $j$ is same as the one from patch $j$ to patch $i$ even if the environment of the two patches are different. Therefore, a symmetric dispersal has no capability to adapt to environmental variations. In the context of (2.1), we may set $c_{i+1i} = c_{ii+1} = \gamma_{i+1/2}$ and rewrite it as

Find that the above is a finite difference scheme for Fick's diffusion law,

In several space dimensions, the equation is written as $u_t = \nabla\cdot(\gamma\nabla u)$. Any constant state is a steady state of (2.2). Remember that a constant steady state is a result of symmetry, not randomness.

Random diffusion. The population in path $i$ may migrate to patch $i+1$ or $i-1$. The dispersal in (2.1) is called random if the two migration rates are identical, i.e., if $c_{i+1i} = c_{i-1i}$. These departing rates may depend on the living condition in patch $i$. In starvation driven dispersal, we assume that the migration rate increases if the environment of the patch becomes harsh. Hence, a spatially heterogeneous random dispersal is automatically asymmetric, i.e., $c_{ii+1}\ne c_{i+1i}$, if the environment of the two patches are different. This asymmetry is an important ingredient as a biological dispersal since it gives the species an ability to adapt to the environment. Let $s_i$ be an indicator that shows how bad the living condition in path $i$ is. If the migration rate to depart path $i$ is decided by this indicator $s_i$, i.e., $c_{i\pm1i} = \gamma(s_i)$ for some increasing function $\gamma(s)$, we may write (2.1) as

This is a finite difference scheme of

which is written as $u_t = \Delta (\gamma(s) u)$ in several space dimensions. A constant distribution is not a steady state anymore. The randomness produces nonconstant steady state if the environment is heterogeneous. The diffusion equation (2.3) is often called Fokker-Planck type diffusion since it contains an advection phenomenon inside. However, asymmetric diffusion is a more specific name for this diffusion since the asymmetry is the key of this diffusion.

Symmetric random diffusion. If the dispersal in (2.1) is random and symmetric at the same time, the dispersal rates are all identical, i.e., $c_{ij} = \gamma_0$ for some constant $\gamma_0 > 0$. Then, (2.1) is written as

This is a finite difference scheme for the heat equation,

which is written as $u_t = \gamma_0\Delta u$ in several space dimensions. In other words, the classical heat equation (2.4) with constant coefficient is a special case with symmetry and randomness. People often incorrectly say that randomness makes the density distribution converges to a constant steady state. It is actually the symmetry of the dispersal that makes a constant distribution a steady state.

3.   Advection related to logistic population dynamics

In this section, we derive advection models when there is only one species and the population dynamics is given by the logistic equation. We will consider two approaches which will respectively give a minimal and a logarithmic type advection terms. We will also find meaningful parameter regimes of these models. The properties of driven dispersal models are compared by testing the structure of steady state solutions and the evolution of dispersal strategies through a competition of two different strategies. The derivation process is a simplification of a starvation driven diffusion. The goal is to find a simplified model and its parameter regime that keeps key features of the original diffusion operator.

Consider the logistic equation with a starvation driven diffusion,

where $m = m(x) > 0$ is the local carrying capacity at $x\in\Omega$. The starvation measure is a dimensionless quantity given by

where the carrying capacity is the maximum amount of population supported by the environment. Note that the migration rate $\gamma$ is a decreasing function with respect to the specific growth rate $\tau$ and an increasing one with respect to the starvation measure $s$. We will abuse notation by denoting $\gamma(m, u, v)$, $\gamma(\tau)$, or $\gamma(s)$ depending on the context.

We take two hypotheses for $\gamma$. First, it is bounded above and below away from zero, i.e., there exist $\underline \gamma$ and $\overline\gamma$ such that

Second, $\gamma$ is an increasing differentiable function of $s$, i.e.,

The starvation measure is the number of population that should share the unit amount for resource. If population increases or resource dwindles, the survival condition worsens and this ratio increases. It is a natural assumption that the migration rate increases in such cases.

The asymmetric diffusion operator is written as

We can find regular diffusion, self-diffusion, and advection from the three terms in (3.5). Eq. (3.1) is written as

We derive three dispersal models based on this equations.

3.1.   Derivation of a logarithmic model

Rewrite (3.6) by

where

Notice that $\gamma(s)$ and $\gamma'(s)s$ have the same order. We obtain a constant chemosensitivity $\chi(s) = {p\over 1+p}$ if $\gamma(s) = s^p$ for $p > 0$. The two hypotheses (3.3) and (3.4) make the equation uniformly parabolic, i.e., the diffusivity $d(s)$ is uniformly bounded away from zero,

They also give a uniform upper bound of the chemosensitivity $\chi(s)$,

The first simplified dispersal model related to the logistic population dynamics is obtained by treating $d(s)$ and $\chi(s)$ as constants, i.e.,

Since the advection term is written as $\chi_0{u\over m}\nabla m = \chi_0{u}\nabla \ln m$, the obtained advection is often called a logarithmic model. The model is often employed in chemotaxis theory. The size of chemosensitivity $\chi_0$ has no upper bound in the chemotaxis theory. However, in the context of starvation driven diffusion, the above derivation shows that chemosensitivity $\chi_0$ is bounded by one.

Remark 1 (Parameter regime for advection). We have obtained a parameter regime for the advection in (3.9), which is $0 < \chi_0 < 1$. This implies that the parameter $\chi_0$ does not exceed one even if an extreme motility function $\gamma$ is chosen as long as it is monotone. Therefore, if a case with $\chi_0 > 1$ is observed, then we can surely say that this advection is not given by the mechanisms of this paper and it could be an evidence for an active gradient sensing.

Remark 2 (Degenerate diffusion). By taking $d(s)$ as a constant, we have obtained a linear diffusion equation (3.9). If we take $\gamma(s) = {1\over 1+p}s^{p}$ without such simplification, we obtain a nonlinear diffusion equation,

Remember that the original equation (3.10) before simplification has a porous medium equation type degenerate nonlinearity.

Remark 3 (Scale invariance). Scale invariance is a property of objects or laws that do not change when scales of length, mass, or other variables, are multiplied by a common factor. Mathematical equations often model such laws and population models derived in this paper have the property. For example, let $u$ be the solution of (3.9) and $m = g(x)$ be the distribution of the carrying capacity. Replaces $g(x)$ with $cg(x)$ for some $c > 0$ by taking a different scale of units. In such a case, the re-scaled function $cu$ is the solution of (3.9).

Under the zero-flux Neumann boundary condition, the steady state of the logarithmic model (3.9) satisfies

where ${\bf n}$ is the outward unit normal vector on the boundary. The reaction term with linear dispersal gives a unique steady state solution (see ^[19,23] for example).^†

^†However, the monotonicity of $\gamma$ is not enough to give the uniqueness of a steady state solution of the nonlinear problem (3.1) (see [18,Remark 3.1]).

In Figure 2, steady states are given with five different advection coefficients. In this simulation we took $d_0 = 1$, $\Omega = (0, \pi)$, and $m(x) = \sin(x)$. In this and following simulations, we have computed the asymptotic limit of corresponding parabolic problems. Since the steady state solution is unique, the longtime asymptotic limit is the steady state. The case with $\chi_0 = 0$ is the steady state without advection. Notice that the diffusion makes the steady state flat and the logistic population dynamics makes the solution have the maximum at $x = \pi/2$. The case with $\chi_0 = 0.5$ gives a better fitness to $m(x) = \sin(x)$. Indeed, we can observe that the steady state solution converges to $m(x)$ as $\chi_0\to1$. If $\chi_0 = 1$, the steady state gives a perfect fitness $\theta(x) = m(x)$. However, the fitness is worsened if the chemosensitivity goes beyond the meaningful regime, $\chi_0 > 1$. In this regime, the population distribution becomes more concentrated at the maximum point than the concentration of the resource itself. This over-fitting case will turn out to be a bad strategy in a competition context. Note that, due to the scale invariance, the relation between the steady state and $m(x)$ is consistent even if we replace $m(x)$ by $cm(x)$ with $c > 0$.

3.2.   Derivation of a minimal model

Minimal model is often used in chemotaxis theory. We will derive a minimal model for a comparison with the logarithmic model. We will see that the minimal model does not keep scaling invariance and lose the meaningful parameter regime. The minimal model is obtained by including ${1\over m}$ in the chemosensitivity $\chi$, i.e.,

This quantity is not dimensionless anymore. It is not constant neither even if $\gamma$ is a power law. If $\tilde\chi(s, m)$ is replaced with a constant, (3.7) gives a minimal model type advection term, i.e.,

Note that there is no bounded parameter regime for $\tilde\chi_0$ for the minimal model. One might see that the logarithmic model comes more naturally from the asymmetric diffusion than the minimal model does.

Remark 4. Let $u$ be the solution of (3.12). If we replace $m(x)$ by $cm(x)$ and $u$ by $cu$ for a constant $c > 0$, then the equation turns into

which is not the original equation and the constant $c$ appears in the advection term. The lack of scale invariance of the minimal model is the reason for the blowup phenomenon in chemotaxis theory.

The steady state of the minimal model (3.12) satisfies

In Figure 3, steady states are given for the same five cases of Figure 2. Since the minimal model does not have the scale invariance, we have tested two resource distribution functions, $m(x) = \sin x$ and $m(x) = 10\sin x$, for comparison. The case with $\chi_0 = 0$ is the same steady state without advection for both cases. The shape of steady states are similar to the logarithmic model case when $\chi_0$ is large. However, for $\chi_0 < 1$, the two cases are quite different.

In conclusion, the logarithmic model (3.9) is a better choice as a biological dispersal model for several reasons. It shows a consistent behavior with respect to the scale of resource distribution. The reasonable quantity assumed to be constant is $\chi(s): = {\gamma'(s)s\over \gamma(s)+\gamma'(s)s}$, not $\tilde\chi(s, m): = {\gamma'(s)s\over \gamma(s)+\gamma'(s)s}\, {1\over m}$. Furthermore, the logarithmic model gives a meaningful parameter regime which is $0 < \chi_0\le1$.

3.3.   Variation comparison between environment and population

It is often believed that the variation of population distribution is smaller than the variation of environment. Let $m$ be the amount of resource and $u$ be the population of a species that consumes the resource. The scaling law with power (or ratio) $p > 0$ is a relation such as

where $c_0 > 0$ is a constant. This scaling power law is often used to model the relation between the resource and population. After a collection of allometric data, the coefficients $c_0$ and $p$ are interpolated. Indeed, cases with $p\cong{3\over4}$ are widely observed. This means that, even if prey population (or food) increases by 10 times more, the predator population increases approximately by $10^{3/4}$ times more (see Hatton et al. ^[8]). Such scaling laws are explained in terms of the ratio of production over population size (see ^[2]). Even if the observed value of the scaling ratio $p$ has some deviation from three quarters due to various noises and other reasons, the scaling ratio is than one almost every time.

After taking the logarithm on (3.13), we obtain

Therefore, if the population variation is smaller than the environment variation, the scaling ratio $p$ satisfies

where the differentiations is with respect to the space variable. The ratio $p$ computed by (3.14) can give us a good indicator about the admissibility of the model equation. However, this local computation is sensitive at the critical points of $m$ and one need to take an example without a critical point or a one that $m$ and $u$ share the same critical points for a cancellation.

In Figure 4(a), the graphs of the scaling ratio $p$ in (3.14) are given for the five cases corresponding to the examples in Figure 2 when $m(x) = 10\sin(x)$ and $u$ is the solution of the logarithmic model (3.9). We can see that $p\le1$ when $\chi_0\le1$. In other words, the meaningful parameter regime for advection term, $\chi_0 < 1$, is equivalent to the scaling condition $p < 1$. These graphs do not change even if we change the scale of the resource distribution, say $m(x) = \sin(x)$. We can conclude that the meaningful parameter regime in this paper gives the same regime that the population variation is smaller than the environment variation.

In Figure 4(b), similar graphs are given when $m(x) = 10\sin(x)$ and $u$ is the solution of the minimal model (3.12). We can see that $p < 1$ only for the case with $\chi_0 = 0$. Furthermore, $p$ has maximum at the critical point of $m$. These graphs depend on the scale of $m$.

3.4.   Evolution of dispersal strategies

One way to see the property of a biological dispersal strategy is to let it compete with other strategies. To see the effect of dispersal strategies more clearly, other factors are usually set identical. We first consider a competition model,

where the super index in $u^i$ denote the $i$-th species. The second species takes the minimal model type advection. This case has been studied by many authors. In particular, a comprehensive picture of the dynamics has been obtained by Lam and Ni ^[19]. For example, if $\chi_2 = 0$ and $d_1 < 1$, then $u^1$ is selected. If $\chi_2 = 0$ and $d_1 > 1$, then $u^2$ is selected. The steady states of this competition model are given in Figure 5 with $d_1 = 1$ and $\chi_2 > 0$. In these simulations, we observe that, if $\chi_2\gg1$ is large, the species with advection agammaegates near the maximum point of $m(x)$ and the other species takes most of the resources as shown in Figures 5(e) and 5(f). One can also observe that the second species $u^2$ is selected and the first one becomes extinct if $\chi_2\lesssim2$ as in Figures 5(a, b, c). However, the two species coexist otherwise.

To obtain a comprehensive picture of extinction and coexistence parameter regimes, the total population of steady states for each $(\chi_2, d_1)\in [0, 4]^2$ are computed numerically and then their contour maps on $(\chi_2, d_1)$ plane are given in Figure 6 using two different resource distributions. The first row is with $m(x) = \sin(x)$ and the second one is with $m(x) = 2\sin(x)$. One can observe the difference between the two rows and conclude that the extinction and coexistence of competing species do not depend only on the parameters $d_1$ and $\chi_2$, but also on the scale of the resource distribution. The figures in the third column are contour maps of the product of two populations, which show the coexistence regime clearly.

Next we consider a competition with the logistic advection model,

where the third species $u^3$ has the logarithmic model type advection. The steady states of this competition system are given in Figure 7 with fixed $d_1 = 1$ and various $\chi_3 > 0$. If $\chi_3\to\infty$, the species with the advection agammaegates toward the maximum point of $m(x)$ and the other species spreads over the domain, which is the same asymptotic behavior as the previous case in Figure 5. It is clearly observed that the meaningful parameter regime $\chi_3 < 1$ gives the criterion for the selection of the logarithmic model type advection. If $\chi_3 > 1$, then the two species coexist.

To obtain a comprehensive picture of the competition dynamics, the total population of steady states are similarly computed for each $(\chi_3, d_1)\in [0, 3]^2$ and their contour maps on the parameter plane are given in Figure 8. We observe from Figure 8(c) that two species coexist if $\chi_3 > 1$. This is the regime beyond the meaningful parameter regime. This indicates the importance of the meaningful parameter regime as a migration strategy of biological organisms. However, even the third species with the logarithmic model type advection becomes extinct if $d_1 > 0$ and $\chi_3 > 0$ are small enough as given in Figure 8(b). It seems that there exists a curve that connects $(\chi_3, d_1) = (1, 0)$ and $(0, 1)$ which separates the region $\chi_3 < 1$ in Figure 8(c) into two subregions. Then, $u^1$ becomes extinct in one side and $u^3$ in the other. The resource distribution $m(x) = \sin(x)$ has been used in the computation of Figure 8. However, for any resource distributions such as $m(x) = k\sin(x)$ with $k > 0$, we obtain the same contour map. The extinction and the coexistence depend on the parameters of the model, but not on the resource distribution.

3.5.   Examples for a different population dynamics

We will be distracted a little bit to do an extra exercise of deriving dispersal model when a different population dynamics is given. Consider a variation of the logistic population dynamics given by

This population dynamics is also often called logistic and have been studied well. However, in this model, $m(x)$ plays the role of growth rate and carrying capacity at the same time. In this example, we will show how to work directly with the specific growth rate,

Rewrite (3.17) as

We will apply the same approach of previous sections. First, we derive an advection model without self-diffusion. To do that we rewrite the equation as

The coefficient ${-\gamma'(\tau)u\over\gamma(\tau)-\gamma'(\tau)u}$ is a dimensionless quantity since $u$ and $\tau$ have the same dimension in the model ($\tau = m-u$). We similarly set $\gamma(\tau)-u\gamma'(\tau) = d_0$ and ${-\gamma'(\tau)u\over\gamma(\tau)-\gamma'(\tau)u} = \chi_0$. Then, we obtain a simplified advection model and a meaningful parameter regime,

One might feel uncomfortable with the advection term since $u$ is not multiplied to it. Remember that this dispersal model is valid only with the population dynamics term in the equation.

In Figure 9, the steady states are given when $d_0 = 1$ and $m(x) = \sin(x)$. Note that the meaningful parameter regime for the model equation (3.19) is $0 < \chi_0\le1$. The steady state with $\chi_0 = 0$ is the case without advection. We can see that the steady states is similar to the case of (3.9) for $\chi_0\le1$ and, as $\chi_0\to1$, they converge to $m(x)$. For $\chi_0 > 1$, the steady states show completely different structure and starts to take negative values. Therefore, the model is meaningless as a biological dispersal model if $\chi_0 > 1$.

Next, we derive an equation with self-diffusion. It seems that the dispersal model with a self-diffusion term is a better choice if the reaction function is given by $u(m(x)-u)$. Rewrite (3.18) as

In this case, the coefficients depend only on $\tau$ and we set $\gamma(\tau) = d_0 > 0$ and $p(\tau) = p_0 > 0$ to be constants. Then, the corresponding simplified equation is

In this model we have self-diffusion and advection terms which share the same coefficient, $p_0u$. If we take $p_0\to\infty$, self-diffusion and the advection terms grow together. This corresponds to the case of (3.19) when $\chi_0\to1$. Note that the appearance of $\nabla(m-u)$ in the equation is not because the organisms sense that the gradient of the difference $u-m$, but because they change the dispersal rate depending on the specific growth rate, which is $m-u$ in the model.

In Figure 9, steady states of (3.20) are given. The meaningful parameter regime for the equation (3.20) is $p_0 > 0$ and there is no upper bound. As $p_0\to\infty$, we can observe that the steady state converges to $m(x)$.

Remark 5. Another widely used population dynamics is $f(u) = m(x)u(1-u)$. The advection and self-diffusion terms can be similarly derived using the specific growth rate $\tau = m(x)(1-u)$ which is left for readers.

4.   Cross-diffusion in Lotka-Volterra competition equations

Dispersal phenomena caused by inter-species interaction are often modeled by cross-diffusion. A well-known cross-diffusion model is the one suggested by Shigesada, Kawasaki and Teramoto, which is often called the SKT model. The starvation driven diffusion also contains cross-diffusion phenomenon and we compare them in this section.

4.1.   Cross-diffusion in SKT model

In the context of the Lotka-Volterra competition system, the SKT model is written as

To save space, we will mostly write the equation for the first species $u^1$ only. The equation for $u^2$ is similarly written. The coefficients are usually normalized by setting $b_{11} = b_{22} = 1$. The SKT model is a random dispersal with a migration rate $\gamma^i = d_i+a_{i1}u^1+a_{i2}u^2$. The coefficients $a_{ij}$'s are given independently from the population dynamics. However, the response of an individual to the environment and other species should be reflected in the coefficients. The key in modeling a cross-diffusion term is in finding meaningful coefficients from the population dynamics.

The diffusion operator of the SKT model is written as

which contains four components. The first three terms are respectively called Fickian diffusion, self-diffusion, and cross-diffusion. The fourth one, $a_{12}u^2\nabla u^1$, may look unusual since it is in the direction of $\nabla u^1$ with a coefficient proportional to the population of the other species.

4.2.   Cross-diffusion in SDD

We derive cross-diffusion terms for a competition system. The coefficients of cross-diffusion terms are decided in terms of coefficients of the reaction term. In the context of Lotka-Volterra competition system in (4.1), the starvation measure for the $i$-th species is given by

This definition is consistent to earlier cases. We consider two species case with $i = 1$ or $2$. To reduce the use of indexes, we denote $u^1, u^2, m^1, s^1, \gamma_1$ by $u, v, m, s, \gamma$, respectively. Then, the starvation driven diffusion for $u$ is written as

This relation is similar to (3.6) except the cross-diffusion term, $\gamma'(s){b_{12}u\over m}\nabla v$.

For simplification, we first rewrite the equation as

where

These two coefficients are the ones corresponding to (3.8). The difference is that $s$ is not ${u\over m}$ anymore. This difference makes the derivation subtle. Note that $\chi$ is still bounded by $0$ and $1$ and we obtain

However, even if $u = 0$, there can be non-zero flux and hence the nonnegativity of a solution might break. This happened in the simplification process and we need to set the flux to be zero when $u = 0$. Hence, we take

Next, we derive a cross-diffusion model with self-diffusion which turns out to be a more natural choice. First, rewrite (4.2) as

where

If we set both $d(s)$ and $p(s)$ be constant, the resulting equation is

where $d_0, p_0 > 0$. Find that the same coefficient $p_0$ appears in the self-diffusion, cross-diffusion, and advection terms. If $u = 0$, both reaction and diffusion terms in (4.5) becomes zero and the solution becomes nonnegative.

Instead of taking constant coefficients, we may keep a nonlinearity that comes from $\gamma$. For example, if we may simply take $\gamma(s) = s$, then (4.2) with a competition term is written as

In this case the diffusion is degenerate.

Remark 6. The driven equations in this section have scaling invariance. For example, if $(u, v)$ is a solution of the problem for a given $m$, its scalar multiplication $(cu, cv)$ with $c > 0$ is a solution when $m$ is replaced by $cm$. The cross diffusion plays a role of advection even if the environment is homogeneous. If we set $m = 1$ in the equations, we obtain equations for the homogeneous case.

4.3.   Competition with cross-diffusion

We investigate the effect of the cross-diffusion in the context of a competition system. In this section, we test the property of cross-diffusion models by numerical computations for the case with $m = 1$ only. In the case the advection term disappears and we can focus on cross-diffusion phenomenon more closely. First, we consider a classical cross-diffusion model,

The space domain is $\Omega = (0, 1)$ and the boundary condition is Neumann as earlier. The coefficients of reaction terms are given by a symmetric matrix:

This system is for a strong competition and the two species are identical in terms of the population dynamics. Hence, the competition is between two different dispersal strategies. The first species $u$ has the linear diffusion only with diffusivity $d_1$. The second species $v$ has the linear diffusion with fixed diffusivity $d = 1$ and a classical cross-diffusion with diffusivity $2\chi v$. We first test the behavior of this classical model. For a strong competition system, the result of competition depends on the choice of the initial value. For a fair comparison we take initial values identically,

In the first row of Figure 10, the competition results are given for three cases of initial values. Since it is a strong competition case, one of the two species goes extinct which is denoted in the figures. Under the initial values satisfying (4.7), the competition result depends on the initial total mass. The first case, $u_0(x) = v_0(x) = {1\over2}x$, is when the initial population size is smaller than the carrying capacity and the population size grows. In this growing stage, the cross-diffusion has positive effect and is selected in a larger parameter regime. The second case, $u_0(x) = v_0(x) = x$, is when the initial population size is same as the carrying capacity. In this case, $v$ is selected in a smaller parameter regime. The last case, $u_0(x) = v_0(x) = 2x$, is when the initial population size is larger than the carrying capacity. In this shrinking stage, the cross-diffusion has negative effect and is selected in an even smaller parameter regime. If two initial values have the same total mass, extinction parameter regimes are similar.^‡

^‡The wiggles along the boundaries are due to the mesh size. The boundaries of the extinction parameter regimes seem to converge to lines as the mesh size decreases.

In the second row of Figure 10, extinction parameter regimes are given when the competition system is replaced by

The cross-diffusion model (4.4) is taken in this system, which has been derived in this paper from the starvation driven diffusion. From the figures in the second row, we can observe that the extinction parameter regimes are relatively invariant under the choice of initial values even if initial total population sizes are different. The parameter regime, in which the species $u$ goes extinct, is bounded by two lines of slope $b_{21}$ which is 2 in the example. The distance of the two boundaries is about $1$ in the direction of $\chi$-axis. Note that the size of meaningful parameter regime of $\chi$ is 1 which is given in (4.3). In other words, the meaningful parameter regime for $\chi$ is exactly the parameter regime that the cross-diffusion is beneficial.

5.   Cross-diffusion in Lotka-Volterra predator-prey equations

Consider Lotka-Volterra predator-prey equations,

where $u$ and $w$ are population densities of prey and predator species, respectively. We put the coefficients $a, b, c$ and $d$ in the model to observe their dependency. We first consider the prey species. The specific growth rate is $\tau_u = a-bw$. Let $s = bw$. Then, $\gamma(s)$ is an increasing function and the prey equation is written as

Let $\gamma(s) = s^p$ for $p > 0$. Then, we have ${\gamma'(s)s\over\gamma} = p$. Now we disregard the common diffusivity by letting $\gamma = d_u$ constant and obtain

Similarly, we may obtain a corresponding equation for predator species $w$;

Notice that the obtained cross-diffusion term is independent of the coefficients $a, b, c, d$, and cross-diffusion models are in logarithmic type.

5.1.   Competition with cross-diffusion

Murray ^[21] showed that the solution of the Lotka-Volterra predator-prey equations converges to a spatially homogeneous solution. Therefore, the effect of cross-diffusion will eventually disappear and testing its property becomes tricky. To investigate the effect of cross-diffusion in a simplest setting we first consider a competition system of two prey species with time independent distribution of predator species. Consider

where the space domain is $\Omega = (0, 2\pi)$. In this system $u$ and $v$ are prey species. The predator species is considered as a given function $w = 1+0.9\sin(x)$. Competition terms are added with carrying capacity $K = 100$ to provide stability. Since the carrying capacity is large enough, the population is controlled by predators. The parameters $p$ and $\chi$ are the sizes of the cross-diffusion terms. Here, the second species $v$ follows the cross-diffusion model derived from the population dynamics. The extinction and coexistence regimes of the two species are given in Figure 11. When $p = 0$, $v$ is selected in the regime of $0 < \chi < 1$, and the two species coexist in the regime $\chi > 1$. As the parameter $p$ increases, these regimes change.

Next we consider a competition system of two predator species with time independent distribution of prey species. Consider

where the space domain is $\Omega = (0, 2\pi)$. In this system $v$ and $w$ are predator species. The prey species is a given function $u = 1+0.9\sin(x)$. This ratio dependent model is considered since the prey should be shared among predators. The ratio ${u\over v+w}$ is the average amount of prey that each predator can take in average. Hence, the population growth of species $v$ is proportional to ${uv\over v+w}$. In Figure 11, the extinction and coexistence regimes of the two species are given. When $\chi = 1$, $v$ is selected for all $p > 0$. The species $v$ becomes extinct in two disjoint regimes when $\chi$ is small or large. The other species $w$ becomes extinct when $\chi < 1$ and $p$ is small or $\chi > 1$ and $p$ is large.

6.   Conclusions

It has been believed for a long time that the random diffusion trivializes initial distribution and the advection or the cross-diffusion phenomena are given by an independent dynamics. However, we have seen in Section 2 that the symmetric diffusion trivializes the initial distribution and random diffusion gives nonconstant steady state if the environment is heterogeneous. Such a random diffusion is given by an asymmetric diffusion (or Fokker-Plank type diffusion) which contains advection and cross-diffusion in a natural way as in (1.1). In this paper we develop a systematic approach to derive dispersal models from a reaction term which reflects the interactions between individuals and responses to the environment.

Many of dispersal models for biological organisms are based on dispersal theories designed for inorganic substances and the parameters are given independently from population dynamics. Because of that, the essence of biological dispersal is often missing. Biological organisms have the ability to adapt to the spatial heterogeneity and temporal fluctuation of environment and biological dispersal should contain such ingredients since the adaptation is the reason for migration. We have seen in the paper that the dynamics to adjust to the environmental variations can be obtained by changing migration rate depending on the harshness of the environment.

In this paper we have developed and tested biological dispersal models based on an asymmetric diffusion equation,

where $\gamma$ is the departing probability of a species and $f$ is the population dynamics. It was assumed that $\gamma = \gamma(\tau)$ is a decreasing function of the specific growth rate $\tau = {f/u}$ or an increasing function of starvation measure $s$ (see (3.2) for example). After a series of simplification processes, this simple assumption lead us to various advection and cross-diffusion models, (1.3)-(1.10), together with meaningful parameter regimes. Numerical simulations show that the extinction and coexistence parameter regimes of are often decided by these meaningful parameter regimes.

Acknowledgments

Authors greatly appreciate valuable comments and suggestions of anonymous reviewers of this paper which significantly improved the initial version of the paper. This research was supported in part by National Research Foundation of Korea (NRF-2017R1A2B2010398).

Conflict of interest

The authors declare that there is no conflict of interest in this paper.