1. Introduction

This paper is inspired by the nano-ecology experiments on the bacterium Escherichia coli by Keymer et al. reported upon in [6] and the subsequent chemotactic reaction-diffusion model developed and analyzed in [3] prompted by the experiments. In the experiment, Keymer et al. fabricated a one-dimensional array of 85 "microhabitat patches (MHP's)" for the E.coli which are linked by corridors. The corridors allow the bacteria to move from one MHP to the next, and are sufficiently more narrow than the MHP's that the MHP's may be viewed as the "nodes" in the overall environment. Each MHP is connected to feeding channels through which a controllable amount of nutrient passes. Specifically, an individual MHP has dimensions $100\mu m$ (length) $\times 100\mu m$ (width) $\times 30\mu m$ (depth). By comparison, a corridor connecting two adjacent MHP's is $50\mu m$ (length) $\times 5\mu m$ (width) $\times 30\mu m$ (depth). Each MHP is separately attached via two feeding channels to the nutrient source. Each channel has $5$ "nanoslits" of dimensions $20\mu m$ (length) $\times 15\mu m$ (width) $\times 200nm$ (depth) through which nutrient and waste but not bacteria may pass, and the device is such that each "nanoslit" may kept open or closed. In [6], the three primary parameters describing the overall bacterial habitat are the "coupling strength" or flow rate through the corridors, the degree of openness of the feeding channels, and the carrying capacity of the individual MHP's.

In [3], the authors develop a reaction-diffusion-advection model based on the experiments reported in [6], while also taking into account that bacteria often aggregate on the basis of "self-attraction mediated by the excretion of chemoattractants". Their derivation leads to a quasilinear parabolic system for the densities of bacteria and nutrient substrate in a one dimensional habitat, in which bacteria self-aggregate and aggregate in response to nutrient substrate abundance.

In this paper, our aim is to model the situation described in [6] by means of a discrete-diffusion (or patch-island) model for the densities of bacteria and nutrient that also can incorporate bacterial self-aggregation and/or bacterial cueing upon nutrient substrate abundance. We believe that this approach allows us to capture the discrete manner in which micro-habitat patches are embedded in the overall "landscape" that is true to the spirit of [6], while also capturing aggregative behavior as highlighted in [3]. Taking account of the linear structure of the fabricated habitat in [6] leads to "nearest neighbor" dispersal between patches, i.e., bacteria from the $i^{th}$ micro-habitat patch are constrained so as to move directly only to patches $i-1$ and $i+1$. Consequently, our patch model takes on a generalized tridiagonal structure. Indeed, the purely diffusive aspect of dispersal in our model by itself would lead to a cooperative tridiagonal model. The asymptotic predictions ([4] and [8]) of such models are particularly clean, predicting that all orbits tend to equilibria, and these facts inform our analysis to a degree. On the other hand, our modeling and analysis are also informed by several other factors. The most interesting scenarios in the model come when some MHP's are net favorable locales while others are not. In such cases, it is reasonable to think in terms of source-sink population dynamics and of keen interest to determine if the rescue effect arising from dispersal from favorable patches to unfavorable ones promotes persistence throughout the landscape. Consequently, movement rates through corridors are key to a mathematical analysis of the model in much the way that diffusion rates would be in a reaction-diffusion setting. Indeed, discrete diffusion models are often viewed as spatially implicit surrogates for reaction-diffusion models and there is a general expectation that results in the reaction-diffusion setting have parallels in the discrete diffusion framework. Here, in particular, we find strong parallels to reaction-diffusion models with indefinite weight functions in reflecting or closed habitats ([1]; [7]). Moreover, our model is also complicated by the aggregative aspects of dispersal that we include. Specifically, when we incorporate aggregative dispersal into the discrete diffusion framework, we obtain a vector field that is Lipschitz continuous but not continuously differentiable. So we can invoke the existence and uniqueness of solutions to initial value problems for first order systems of ordinary differential equations to guarantee that we have a well-defined dynamical system. However, we cannot directly apply the Hartman-Grobman Theorem to analyze stability of equilibria. We circumvent this obstacle in our asymptotic analysis by drawing upon results from persistence theory ([2]), specifically the Acyclity Theorem and practical persistence estimates.

Two further aspects of our model merit comment at this point. First of all, our patch model has the feature that at equilibrium, if the bacteria density is zero in some patch, it is zero in all patches. That such is the case is an inherent property of the modeling framework. Indeed, the discrete diffusion component of dispersal leads to an "irreducible" coupling of the patches unless some corridor is completely closed, in which case we have two entirely separated landscapes. This feature is analogous to the maximum principle in a reaction-diffusion context. However, one should note that bacteria density in unfavorable patches "far away" from a source patch would be expected to be arbitrarily low and we will illustrate this point via numerical calculations.

The second feature of our model that we should note concerns the dynamics within a single patch in the absence of dispersal between patches. Like [6] and [3], we posit logistic dynamics for the bacteria. The traditional "r-K" form of the logistic model, used in both [6] and [3] has

$\frac{dx}{dt}=rx\left(1-\frac{x}{K}\right),$

(1)

where $x$ is population density, r is the intrinsic growth rate and $K$ is the carrying capacity of the habitat. In [6] and [3] and in our model, the bacterial density is linked to the nutrient density through the dependence of $r$ on nutrient density. In unfavorable patches $r$ will be negative or zero. In such case, there is no mechanism for growth and one should expect that the bacterial density tends toward zero. However, if $r<0$ and $x>K$, $dx/dt>0$. For this reason, we will write the logistic equation in the form

$\frac{dx}{dt}=rx-\alpha x^2,$

(2)

where $\alpha>0$. When $r>0$, (2) is equivalent to (1) with $K=r/\alpha$. When $r\leq 0$, $x(t)$ tends to zero as $t\rightarrow \infty$ for all nonnegative $x(0)$. This alteration has a significant ramification for the coupled bacteria-nutrient interaction within a single patch. Namely, there are only two equilibrium configurations possible for the bacteria-nutrient system as opposed to the three in [6]. The first possibility is that the nutrient influx rate is too low relative to the component $d$ of $r$ that represents the natural mortality of the bacteria and the bacteria population vanishes and the nutrient is at its input value. The other is that nutrient influx is high enough relative to $d$ for the bacteria to survive and the bacteria density is at the $r/\alpha$ value determined by the balance of resource availability and natural mortality. The third equilibrium in [6] results from having a nutrient level in which $r=0$. In this event, damped oscillatory behavior leading to this equilibrium is observed in [6]. By employing (2) in place of (1) we do not see such behavior in the predictions of our model.

The remainder of the paper is structured as follows. We discuss within patch dynamics in section 2, leading to the development of the multi-patch model in section 3. We present our mathematical analysis of the model in sections 4 and 5. In section 4, we consider general discrete diffusion rates, whereas in section 5, we assume diffusion rates are equal. We discuss the results of our numerical investigations in section 6. We end by drawing some biological conclusions in section 7.

---

2. The single patch model

As noted in the introduction, the single patch dynamics of the system in [6] and [3] are described by

$\begin{aligned} \frac{dp}{dt} & =(\mu s-d)p \left(1-\frac{p}{K}\right)\\ \frac{ds}{dt} & =\lambda(1-s)-\varepsilon\mu s p \end{aligned}$

(3)

where $p$ is the density of the bacteria, $s$ is the density of the substrate that the bacteria feed upon, $\mu$ is the bacterial growth rate under the maximal sustainable substrate level and $d$ is bacterial mortality. $K$ represents the carrying capacity of the habitat, $\lambda$ is the flow rate through the feeding corridor, and $\varepsilon$ is a conversion factor from the consumption of substrate to the birth of bacteria. In this formulation, $\varepsilon>1$. As discussed in the introduction, we modify (3) to (4) below to allow for the possibility that $\mu s-d$ can be negative. In so doing, $K$ will no longer represent carrying capacity. The other parameters retain the same meanings as in (3). We make one further adjustment that will be significant once we are in the multi-patch setting. Namely, as scaled by $1-s$, the maximum sustainable density for the substrate is one. As described in [6], the degree of openness of the feeding channels is adjustable. As a result, we want to allow a maximal sustainable density below one in that case. We use the parameter $\beta \in (0,1]$ for this purpose. Consequently, our single patch model is expressed as

$\begin{aligned} \frac{dp}{dt} & = (\mu s-d)p-\frac{p^2}{K}\\ \frac{ds}{dt} & = \lambda(\beta-s)-\varepsilon\mu sp. \end{aligned}$

(4)

To analyze (4), observe first that (4) has either one or two equilibria. To this end, note that the substrate density $s$ is a lower solution to

$\frac{{dz}}{{dt}} = \lambda (\beta - z),$

which converges over time to $\beta$. If $\mu\beta-d<0$, then the bacterial density is a lower solution to

$\frac{{dw}}{{dt}} = (\mu \beta - d)w,$

for $t\gg 0$ and consequently $(p, s)\rightarrow (0, \beta)$ as $t\rightarrow \infty$. A modified argument holds when $\mu\beta-d=0$.

So now suppose $\mu\beta-d>0$. An equilibrium $(p^*, s^*)$ with $p^*\neq 0$ satisfies

${p^*} = K(\mu {s^*} - d),$

and

${p^*} = \frac{{\lambda (\beta - {s^*})}}{{\varepsilon \mu {s^*}}},$

so that $s^*$ must satisfy

$\begin{aligned} \varepsilon\mu^2 K (s^*)^2+(\lambda-\varepsilon\mu Kd)s^*-\lambda\beta=0 \end{aligned}$

(5)

so that

${s^*} = \frac{{(\varepsilon \mu Kd - \lambda ) + \sqrt {{{(\varepsilon \mu Kd - \lambda )}^2} + 4\lambda \beta \varepsilon {\mu ^2}K} }}{{2\varepsilon {\mu ^2}K}}$

is the only positive root of (5). Setting $f(s)=\varepsilon\mu^2 K s^2+(\lambda-\varepsilon\mu Kd)s-\lambda\beta$, we have $f(0)=-\lambda\beta$ while $f(\beta)=\varepsilon\mu K\beta(\mu\beta-d)>0$ so that $s^*<\beta$. Additionally, $f(d/\mu)=\lambda(d/\mu-\beta)<0$ so that $d/\mu<s^*<\beta$ and $0<p^*<K(\mu\beta-d)$. In particular, there is a unique componentwise positive equilibrium to (4).

It is easy to check that the Jacobian matrix for an equilibrium to (4) is

$\left[ {\begin{array}{*{20}{c}} {\mu s - d - \frac{{2p}}{K}} & {\mu p}\\ { - \varepsilon \mu s} & { - \lambda - \varepsilon \mu p} \end{array}} \right]$

which is

$\left[ {\begin{array}{*{20}{c}} {\mu \beta - d} & 0\\ { - \varepsilon \mu \beta } & { - \lambda } \end{array}} \right]$

when $(p,s)=(0,\beta)$. Consequently $(0,\beta)$ is a saddle point unstable in the $p$ direction. Since

$\begin{array}{l} \;\;\frac{{ds}}{{dt}} \le \lambda (\beta - s)\\ \frac{{dp}}{{dt}} \le (\mu \beta - d)p - \frac{{{p^2}}}{K} \end{array}$

for $t$ sufficiently large, orbits to (4) with $p(0)>0$ and $s(0)>0$ are drawn into the positively invariant rectangle

$D = \{ (p,s)|0 \le p \le K(\mu \beta - d),0 \le s \le \beta \} .$

In particular, $(p^*,s^*)$ lies within the rectangle. At $(p^*,s^*)$ the Jacobian matrix simplifies to

$\left[ {\begin{array}{*{20}{c}} { - \frac{{{p^*}}}{K}} & {\mu {p^*}}\\ { - \varepsilon \mu {s^*}} & { - \frac{{\lambda \beta }}{{{s^*}}}} \end{array}} \right]$

Thus the determinant of the Jacobian matrix at $(p^*,s^*)$ is positive and the trace is negative, so that $(p^*,s^*)$ is locally asymptotically stable. Setting

$B = \frac{1}{{ps}},F = (\mu s - d)p - \frac{{{p^2}}}{K},{\rm{and}}\;G = \lambda (\beta - s) - \varepsilon \mu ps$

one may calculate that

$\frac{{\partial (BF)}}{{\partial p}} + \frac{{\partial (BG)}}{{\partial s}} < 0$

in the interior of $D$. Consequently, we have by Dulac's Criterion and Poincar$\acute{e}$-Bendixson theorem that $(p^*,s^*)$ is in fact globally asymptotically stable relative to orbits with $p(0)>0, s(0)>0$.

---

3. The multi-patch model

Here we take a discrete-diffusion (patch island) approach to model dispersal between adjacent patches, which is linear, augmented by a discrete version of nonlinear chemotactic aggregation, where the advection is biased toward higher conspecific density or higher substrate concentration. Both effects are included in the reaction-advection-diffusion model in [3] which mirrors discussion in [6]. The linear level dispersal is standard in models of this type. To model chemotactic self-aggregation, we will posit a tendency for the bacteria to go from, say, patch $i$ to patch $(i+1)$ whenever the bacterial density is higher in patch $(i+1)$ than in patch $i$. Moreover, as the simplest quantitative representation, we will assume this tendency is proportional to the difference between the densities in the two patches. As such, a term of the form

$\gamma {p_i}\max \{ {p_{i + 1}} - {p_i},0\}$

is subtracted from $dp_{i}/dt$ and added to $dp_{i+1}/dt$, where $p_{i}$ is the bacterial density in patch $i$. The corresponding term in modeling chemotactic advection toward higher substrate concentration is

$\nu {p_i}\max \{ {s_{i + 1}} - {s_i},0\}$

where $s_{i}$ is substrate concentration in patch $i$. Once bacteria move from one patch to another, they are identified as residents of the new patch, so that local population dynamics in patch $i$ only involve $p_{i}$ and $s_{i}$. Our analysis will consider the case of an arbitrary number of patches (MHP's). The full model takes the form

$\begin{aligned} \frac{dp_{1}}{dt} & = -D_{21}p_{1}+D_{12}p_{2}-\gamma p_{1} \max\{p_{2}-p_{1},0\}+\gamma p_{2} \max\{p_{1}-p_{2},0\}\\ & \, -\nu p_{1} \max\{s_{2}-s_{1},0\}+\nu p_{2} \max\{s_{1}-s_{2},0\}+(\mu s_{1}-d)p_{1}-\alpha p_{1}^2\\ \frac{ds_{1}}{dt} & = \lambda (\beta_{1}-s_{1})-\varepsilon\mu s_{1}p_{1}\\ \vdots\\ \frac{dp_{i}}{dt} & = D_{i,i-1}p_{i-1}+D_{i,i+1}p_{i+1}-D_{i-1,i}p_{i}-D_{i+1,i}p_{i}\\ & \, +\gamma p_{i-1} \max\{p_{i}-p_{i-1},0\}-\gamma p_{i} \max\{p_{i-1}-p_{i},0\}\\ & \, +\gamma p_{i+1} \max\{p_{i}-p_{i+1},0\}-\gamma p_{i} \max\{p_{i+1}-p_{i},0\}\\ & \, +\nu p_{i-1} \max\{s_{i}-s_{i-1},0\}-\nu p_{i} \max\{s_{i-1}-s_{i},0\}\\ & \, +\nu p_{i+1} \max\{s_{i}-s_{i+1},0\}-\nu p_{i} \max\{s_{i+1}-s_{i},0\}+(\mu s_{i}-d)p_{i}-\alpha p_{i}^2\\ \frac{ds_{i}}{dt} & = \lambda (\beta_{i}-s_{i})-\varepsilon\mu s_{i}p_{i}\\ \vdots\\ \frac{dp_{n}}{dt} & = -D_{n-1,n}p_{n}+D_{n,n-1}p_{n-1}-\gamma p_{n} \max\{p_{n-1}-p_{n},0\}\\ & \, +\gamma p_{n-1} \max\{p_{n}-p_{n-1},0\}-\nu p_{n} \max\{s_{n-1}-s_{n},0\}\\ & \, +\nu p_{n-1} \max\{s_{n}-s_{n-1},0\}+(\mu s_{n}-d)p_{n}-\alpha p_{n}^2\\ \frac{ds_{n}}{dt} & = \lambda (\beta_{n}-s_{n})-\varepsilon\mu s_{n}p_{n} \end{aligned}$

(6)

where $i=2,\ldots,n-1$. Here $D_{ij}$ denotes the rate of diffusion from patch $j$ into patch $i$ (note that $j=i-1$ or $j=i+1$), $\gamma$ denotes the strength of chemotactic bacterial self-aggregation, and $\nu$ the strength of bacterial cueing upon relative substrate concentrations. The single patch parameters $\mu$, $d$, $\lambda$ and $\varepsilon$ have the same meanings as in the preceding section, as does also $\beta_{i}$, except $\beta_{i}$ is allowed to vary from patch to patch. As noted in the introduction, this feature reflects the design of the experimental device in [6]; namely, it is possible to control the maximal available substrate concentration in each patch. Finally, we find it convenient to replace $1/K$ by $\alpha$ in the bacterial logistic dynamics. The parameters $\gamma$, $\nu$, $\mu$, $d$, $\varepsilon$ and $\alpha$ reflect bacterial traits which can be reasonably considered as independent of patch. Similar reasoning applies to $\lambda$ relative to the substrate.

In the analysis that follows, we will set $\nu=0$ so as to focus primarily on bacterial self-aggregation. Our motivation here is several-fold. First of all, when bacteria cue upon conspecific density, it is not unreasonable to think of them indirectly cueing upon relative substrate concentration. Second, a term such as

$\nu {p_i}\max \{ {s_{i + 1}} - {s_i},0\}$

is bounded by $\nu\beta_{i+1} p_{i}$. Hence it is relatively insignificant as a dispersal term compared to the linear movement terms if $\nu$ is small relative to the $D_{ij}$ independent of bacteria density. Such is not the case with the bacterial self-aggregation. Thirdly, our analysis turns to a large extent on the observation that the system with either form of chemotactic aggregation or both is irreducible in the bacterial density. By this statement, we mean that if $p_{i}\equiv 0$ for some $i$, it is identically zero for all $i$. Such follows from the tridiagonal structure and the positivity of diffusion coefficients. This feature can be regarded as akin to a maximum principle in a reaction-diffusion setting. Consequently, we will focus only on the case with bacterial self-aggregation which we believe to be interesting in and of itself.

---

4. Analysis Ⅰ: General diffusion rates

We will give persistence results for the n-patch analogue of (6) with $\nu=0$ under the assumption that the diffusion rates $D_{ij}$ ($j=i-1$ or $i+1$) are positive. Here, it is useful first to present the special case when $n=2$ separately. To this end, we consider the model

$\begin{aligned} \frac{dp_{1}}{dt}= & (-D_{21}p_{1}+D_{12}p_{2})-\gamma p_{1}\max\{p_{2}-p_{1},0\}\\ & +\gamma p_{2}\max\{p_{1}-p_{2},0\}+(\mu s_{1}-d)p_{1}-\alpha p_{1}^2\\ \frac{ds_{1}}{dt}= & \lambda (\beta_{1}-s_{1})-\varepsilon\mu s_{1}p_{1}\\ \frac{dp_{2}}{dt}= & (D_{21}p_{1}-D_{12}p_{2})-\gamma p_{2} \max\{p_{1}-p_{2},0\}\\ & +\gamma p_{1} \max\{p_{2}-p_{1},0\}+(\mu s_{2}-d)p_{2}-\alpha p_{2}^2\\ \frac{ds_{2}}{dt}= & \lambda (\beta_{2}-s_{2})-\varepsilon\mu s_{2}p_{2}\\ \end{aligned}$

(7)

Here we are interested in the solution flow for (7) on the set $X=\{(p_{1},s_{1},p_{2},s_{2}): p_{i}\geq0,0\leq s_{i}\leq\beta_{i},i=1,2\}$. Ultimately, our analysis of (7) on $X$ relies on Theorem 4.5 of [9]. Thieme's result is an extension of the celebrated Acyclicity Theorem of persistence theory due to [5]. A key hypothesis in Theorem 4.5 of [9] is that the solution flow for (7) is forward invariant on the set

${X_1} = \{ ({p_1},{s_1},{p_2},{s_2}):{p_i} > 0,0 < {s_i} < {\beta _i},i = 1,2\} ,$

so we start our discussion with a verification of this fact.

Proposition 1. The solution flow for (7) is forward invariant on $X_{1}$.

Proof. Suppose that $(p_{1}(t),s_{1}(t),p_{2}(t),s_{2}(t))$ is a solution of (7) with $p_{i}(0)>0$ and $0<s_{i}(0)<\beta_{i}$ for $i=1,2$. Suppose that this solution does not remain in $X_{1}$ for all $t>0$. By continuity, there is a smallest $t_{1}>0$ so that $(p_{1}(t),s_{1}(t),p_{2}(t),s_{2}(t))\in X_{1}$ for $0\leq t<t_{1}$ with $(p_{1}(t_{1}),s_{1}(t_{1}),p_{2}(t_{1}),s_{2}(t_{1}))\in \partial X_{1}$. In such case, either $s_{i}(t_{1})=\beta_{i}$ or $0$ for at least one $i$ or $p_{i}(t_{1})=0$ for at least one $i$.

So consider

$\frac{{d{s_i}}}{{dt}} = \lambda ({\beta _i} - {s_i}) - \varepsilon \mu {s_i}{p_i}$

Since $s_{i}\geq0$ and $p_{i}\geq0$ on $[0,t_{1}]$, $s_{i}$ is a subsolution of the problem

$\begin{aligned} \frac{dz}{dt} & =\lambda(\beta_{i}-z)\\ z(0) & =s_{i}(0) \end{aligned}$

on $[0,t_{1}]$. Hence $s_{i}(t_{1})<\beta_{i}$. On the other hand, if we let $M_{i}=\max\{p_{i}(t):t\in[0,t_{1}]\}$, $s_{i}$ is a supersolution of

$\begin{aligned} \frac{dw}{dt} & =\lambda\left(\beta_{i}-\left(1+\frac{\varepsilon\mu M_{i}}{\lambda}\right)w\right)\\ w(0) & =s_{i}(0) \end{aligned}$

which guarantees that $s_{i}(t_{1})>0$. So now consider the equation for $p_{1}$. It is easy to see that $p_{1}$ is a supersolution of

$\begin{aligned} \frac{dy}{dt} & =(\mu s_{1}-d-D_{21}-\gamma \max\{p_{2}-p_{1},0\})y-\alpha y^2\\ y(0) & =p_{1}(0) \end{aligned}$

Thus $p_{1}(t_{1})>0$. An analogous argument showing that $p_{2}(t_{1})>0$ leads us to conclude that $(p_{1}(t_{1}),s_{1}(t_{1}),p_{2}(t_{1}),s_{2}(t_{1}))$ is not in $\partial X_{1}$, a contradiction. So the solution flow is forward invariant on $X_{1}$.

Observe, for example, that if $s_{i}(0)=\beta_{i}$ and $p_{i}(0)>0$, then $s_{i}(t)<\beta_{i}$ for all $t>0$. Consequently, the solution flow of (7) is not forward invariant on $\partial X_{1}$. The Acyclicity Theorem of [5] would require such. It is one of the key insights of [9] to remove this requirement.

Next we show that the solution flow of (7) is asymptotically bounded or point dissipative. To this end, observe that

$\begin{aligned} \frac{d}{dt}(p_{1}+p_{2}) & =(\mu s_{1}-d)p_{1}-\alpha p_{1}^2+(\mu s_{2}-d)p_{2}-\alpha p_{2}^2\\ & \leq \rho(p_{1}+p_{2})-\alpha(p_{1}^2+p_{2}^2) \end{aligned}$

where $\rho=\max\{\mu\beta_{1}-d,\mu\beta_{2}-d,1\}$. Now $(p_{1}+p_{2})^2\leq 2(p_{1}^2+p_{2}^2)$, which implies that

$\rho(p_{1}+p_{2})-\alpha(p_{1}^2+p_{2}^2)\leq \rho(p_{1}+p_{2})-\frac{\alpha}{2}(p_{1}+p_{2})^2.$

Consequently, given any $\sigma>0$, there is a $T=T(p_{1}(0),p_{2}(0))$ such that

$p_{1}(t)+p_{2}(t)\leq \frac{2\rho}{\alpha}+\sigma$

for $t\geq T$. In light of Proposition 1, this establishes

Proposition 2. Solutions of (7) are asymptotically bounded in $X$.

We observe that Propositions 1 and 2 extend in an analogous manner to the n-patch analogue of (6) with $\nu=0$. We shall use this fact in the sequel as needed without further argument.

Verifying permanence or uniform persistence in (7) via Theorem 4.5 of [9] requires an understanding of the flow of (7) restricted to the boundary of $X_{1}$. Here effectively the boundary of $X_{1}$ is having $s_{i}\equiv\beta_{i}$ or $p_{i}\equiv0$ for $i=1$ or $i=2$. It is immediate that $s_{i}\equiv\beta_{i}$ for $i=1$ or $2$ implies that $p_{i}\equiv0$. The tridiagonal structure of the model at the linear level forces $p_{j}\equiv0$ where $j\neq i$. Having $p_{j}\equiv0$ in turn implies that $s_{j}\rightarrow\beta_{j}$ as $t\rightarrow\infty$. Consequently, the only element of $w(\partial X_{1})$ is $\{(0,\beta_{1},0,\beta_{2})\}$ and the acyclicity requirement in Theorem 4.5 of [9] is automatic (again due to the tridiagonal structure, which is analogous to having a maximum principle in a reaction-diffusion model). So it becomes the case that (7) is permanent or uniformly persistent if and only if the stable manifold of $\{(0,\beta_{1},0,\beta_{2})\}$, denoted $W^{s}(\{(0,\beta_{1},0,\beta_{2})\})\cap X_{1}=\emptyset$. So we have:

Theorem 4.1. The system (7) is permanent (uniformly persistent) if and only if

$W^{s}(\{(0,\beta_{1},0,\beta_{2})\})\cap X_{1}=\emptyset.$

We now examine when Theorem 4.1 holds. We begin with some simple observations. Note from (7) that

$\frac{d(p_{1}+p_{2})}{dt}=(\mu s_{1}-d)p_{1}+(\mu s_{2}-d)p_{2}-\alpha (p_{1}^2+p_{2}^2)$

(8)

If $\mu\beta_{1}-d<0$ and $\mu\beta_{2}-d<0$, we have by (8) that $p_{1}+p_{2}$ is a subsolution of the initial value problem

$\begin{array}{l} \frac{{dz}}{{dt}} = z\\ z(0) = {p_1}(0) + {p_2}(0) \end{array}$

where $\underline{c}=\max\{\mu\beta_{1}-d,\mu\beta_{2}-d\}<0$. Consequently, for any initial data with $p_{i}(0)\geq0$, we get that $p_{1}\rightarrow0$ and $p_{2}\rightarrow0$ as $t\rightarrow\infty$, which in turn implies that $W^{s}(\{(0,\beta_{1},0,\beta_{2})\})\cap X_{1}\neq\emptyset$. On the other hand, suppose that $\mu\beta_{1}-d>0$ and $\mu\beta_{2}-d>0$ and that $(p_{1}(t),s_{1}(t),p_{2}(t),s_{2}(t))\rightarrow(0,\beta_{1},0,\beta_{2})$ for some initial data $(p_{1}(0),s_{1}(0),p_{2}(0),s_{2}(0))$ with $p_{i}(0)\geq0$ and $p_{1}(0)+p_{2}(0)>0$. Then $p_{i}(t)>0$ and $0<s_{i}(t)<\beta_{i}$ for all $t>0$ and for $t$ sufficiently large, say $t\geq T_{1}$, $p_{1}+p_{2}$ is a supersolution of the initial value problem

$\begin{array}{l} \frac{{dz}}{{dt}} = cz - \alpha {z^2}\\ z({T_1}) = {p_1}({T_1}) + {p_2}({T_1}), \end{array}$

where $c$ is positive and less than $\min\{\mu\beta_{1}-d,\mu\beta_{2}-d\}$.

Consequently, given any $\sigma\in(0,c/\alpha)$, there is a $T_{2}>T_{1}$ so that $p_{1}(t)+p_{2}(t)>c/\alpha-\sigma$ for all $t\geq T_{2}$. In particular, it is not possible for a solution of (7) starting in $X_{1}$ to converge to $\{(0,\beta_{1},0,\beta_{2})\}$. So if $\mu\beta_{1}-d>0$ and $\mu\beta_{2}-d>0$, $W^{s}(\{(0,\beta_{1},0,\beta_{2})\})\cap X_{1}=\emptyset$. The upshot is as follows. If both microhabitat patches are unfavorable, the bacteria goes extinct in both patches over time, regardless of its dispersal pattern. On the other hand, if both patches are favorable, the bacterial species persists in both, again regardless of its dispersal pattern. Consequently, the only case in which dispersal behavior can have an impact on the asymptotic predictions of the model is when we have so-called source sink dynamics, where in the growth rate $\mu\beta_{i}-d$ is positive in one patch (meaning a favorable or source habitat) but negative in the other (meaning an unfavorable or sink habitat).

So for the remainder of this section, we assume that $\mu\beta_{1}-d>0$ while $\mu\beta_{2}-d<0$. Assume initially that there is no bacterial self-aggregation; i.e., $\gamma=0$ in (7). In this case, the Hartman-Grobman Theorem is applicable to enable us to decide whether or not $W^{s}(\{(0,\beta_{1},0,\beta_{2})\})\cap X_{1}=\emptyset$.

The Jacobi matrix at $(0,\beta_{1},0,\beta_{2})$, $J(0,\beta_{1},0,\beta_{2})$ is given by

$J(0,\beta_{1},0,\beta_{2})= \begin{bmatrix} -D_{21}+\mu\beta_{1}-d & 0 & D_{12} & 0\\ -\varepsilon\mu\beta_{1} & -\lambda & 0 & 0\\ D_{21} & 0 & -D_{12}+\mu\beta_{2}-d & 0\\ 0 & 0 & -\varepsilon\mu\beta_{2} & -\lambda \end{bmatrix}$

(9)

We will use $\sigma$ to denote the eigenvalues of $J(0,\beta_{1},0,\beta_{2})$ in (9). It is not difficult to see that $\sigma=-\lambda$ is a double eigenvalue (reflecting the restorative tendency of the substrate when bacterial abundance is low) and the other eigenvalues of $J(0,\beta_{1},0,\beta_{2})$ are the eigenvalues of the matrix $A$ given by

$A= \begin{bmatrix} -D_{21}+(\mu\beta_{1}-d) & D_{12}\\ D_{21} & -D_{12}+(\mu\beta_{2}-d) \end{bmatrix}.$

(10)

So the other two eigenvalues are given by

$\sigma=\frac{Tr\pm\sqrt{(Tr)^2-4\Delta}}{2}$

where $Tr=\mu(\beta_{1}+\beta_{2})-2d-(D_{12}+D_{21})$ and $\Delta=(\mu\beta_{1}-d)(\mu\beta_{2}-d)-D_{12}(\mu\beta_{1}-d)-D_{21}(\mu\beta_{2}-d)$ are the trace and determinant of $A$. When $D_{12}=D_{21}=0$, there is no dispersal between microhabitat patches and we have persistence of the bacteria in patch 1 and extinction in patch 2. Notice that when $D_{12}=0=D_{21}$, $\Delta<0$ while $Tr$ could be of either sign. When $D_{12}$ and $D_{21}$ are positive and small, $\Delta$ remains negative. Consequently, $J(0,\beta_{1},0,\beta_{2})$ admits one positive eigenvalue and we have that $W^{s}(\{(0,\beta_{1},0,\beta_{2})\})\cap X_{1}=\emptyset$, and the model predicts persistence of the bacteria in both microhabitat patches. Here the subsidy from patch one rescues the bacterial population in patch two, an example of source-sink dynamics at work. If $Tr<0$ when $D_{12}=0=D_{21}$, it remains negative for all positive values of $D_{12}$ and $D_{21}$. If $Tr>0$ when $D_{12}=0=D_{21}$, it will eventually become negative as $D_{21}$ increases. By re-writing $\Delta$ and $Tr$ as

$\begin{aligned} \Delta & = (\mu\beta_{1}-d-D_{21})(\mu\beta_{2}-d)-D_{12}(\mu\beta_{1}-d)\\ Tr & = (\mu\beta_{1}-d-D_{21})+(\mu\beta_{2}-d-D_{12}) \end{aligned}$

it is easy to see that $\Delta$ remains negative as $D_{21}$ increases until after the point at which $Tr$ becomes negative as $D_{21}$ increases. It is also easy to see from the above that $\Delta$ will become positive once $D_{21}$ is large enough. When this happens both eigenvalues of $A$ have negative real parts and we conclude that $W^{s}(\{(0,\beta_{1},0,\beta_{2})\})\cap X_{1}\neq\emptyset$. Consequently, the bacteria from the first patch over disperse to the second patch, effectively turning the first patch into a sink. The bacteria go extinct in both patches if the density in patch $1$ becomes too diminished.

So suppose now that $D_{12}>0$, $D_{21}>0$ and $\gamma>0$, so that we take bacterial self-aggregation into account. Theorem 4.1 remains valid. Suppose that $D_{21}$ is small enough so that $\mu\beta_{1}-d-D_{21}>0$. Let $r_{1}$ and $r_{2}$ be positive numbers so that

$\mu(\beta_{1}-r_{1})-d-D_{21}-r_{2}>0$

and suppose there is a solution of (7) with $p_{i}(0)>0$, $s_{i}(0)\in(0,\beta_{i}),i=1,2$ so that $(p_{1}(t),s_{1}(t),p_{2}(t),$ $s_{2}(t))\rightarrow(0,\beta_{1},0,\beta_{2})$ as $t\rightarrow\infty$. Choose $T>0$ so that for $t\geq T$,

$\begin{aligned} p_{1}(t) & < w<\frac{\mu(\beta_{1}-r_{1})-d-D_{21}-r_{2}}{\alpha}\\ s_{1}(t) & > \beta_{1}-r_{1}\\ \end{aligned}$

and

$p_{2}(t) < \frac{r_{2}}{\gamma}.$

Then

$\begin{aligned} \frac{dp_{1}}{dt}= & (-D_{21}p_{1}+D_{12}p_{2})-\gamma p_{1}\max\{p_{2}-p_{1},0\}+\gamma p_{2}\max\{p_{1}-p_{2},0\}\\ & +(\mu s_{1}-d)p_{1}-\alpha p_{1}^2\\ \geq & -D_{21}p_{1}-\gamma p_{1} \max\{p_{2}-p_{1},0\}+(\mu s_{1}-d)p_{1}-\alpha p_{1}^2. \end{aligned}$

Now $\max\{p_{2}-p_{1},0\}\leq p_{2}$. So $-\gamma p_{1} \max\{p_{2}-p_{1},0\}\geq -\gamma p_{1}p_{2}$. Since $p_{2}<r_{2}/\gamma$, we get that

$\begin{aligned} -\gamma p_{1}\max\{p_{2}-p_{1},0\}\geq -r_{2}p_{1}. \end{aligned}$

Hence we obtain that for $t\geq T$,

$\frac{{d{p_1}}}{{dt}} \ge (\mu ({\beta _1} - {r_1}) - d - {D_{21}} - {r_2}){p_1} - \alpha p_1^2.$

Consequently, $p_{1}(t)\geq \rho(t)$, where $\rho(t)$ is the solution of

$\begin{aligned} \frac{d\rho}{dt}=(\mu(\beta_{1}-r_{1})-d-D_{21}-r_{2})\rho-\alpha\rho^2 \end{aligned}$

on $(T,\infty)$ with $\rho(T)=p_{1}(T)$. Since $\rho\rightarrow\frac{\mu(\beta_{1}-r_{1})-d-D_{21}-r_{2}}{\alpha}$ as $t\rightarrow\infty$, we get that $p_{1}(t)>w$ for $t\gg T$, a contradiction. So there can be no such orbit and thus

${W^s}(\{ (0,{\beta _1},0,{\beta _2})\} ) \cap {X_1} = \emptyset .$

So (7) remains permanent when $\gamma>0$, $D_{12}>0$ and $D_{21}$ is positive and small.

When $\gamma=0$ and $D_{12}>0$ is fixed, (7) loses permanence when $D_{21}$ becomes large enough so that the determinant $\Delta$ of $A$ in (10) becomes positive. For $\gamma>0$ and $D_{12}>0$ and fixed, we establish that (7) is not permanent for $D_{21}$ sufficiently large via the following result.

Proposition 3. For $\gamma>0$ and $D_{12}>0$ fixed, (7) does not admit a componentwise positive equilibrium for $D_{21}$ sufficiently large.

Proof. Suppose (7) admits an equilibrium with $p_{i}>0$ and $s_{i}\in(0,\beta_{i})$ for $i=1,2$. Then we have that

$D_{21}p_{1}+\gamma p_{1}\max\{p_{2}-p_{1},0\}-(\mu s_{1}-d)p_{1}+\alpha p_{1}^2=(D_{12}+\gamma \max\{p_{1}-p_{2},0\})p_{2}$

which is equivalent to

$\begin{aligned} D_{21}+\gamma \max\{p_{2}-p_{1},0\}-(\mu s_{1}-d)+\alpha p_{1}\\ =(D_{12}+\gamma \max\{p_{1}-p_{2},0\})\left(\frac{p_{2}}{p_{1}}\right) \end{aligned}$

(11)

For the time dependent problem, we have

$\frac{d}{dt}(p_{1}+p_{2})\leq (\mu\beta_{1}-d)(p_{1}+p_{2})-\frac{\alpha}{2}(p_{1}+p_{2})^2$

so that $p_{1}$ and $p_{2}$ are ultimately bounded above by any number greater than $[2(\mu\beta_{1}-d)/\alpha]$, independent of the value of $D_{21}$. Therefore, componentwise positive equilibria are bounded independent of $D_{21}$.

One easily observes that the left hand side of (11) tends to $+\infty$ as $D_{21}\rightarrow +\infty$. So if we have componentwise positive equilibria to (7) for arbitrarily large values of $D_{21}$, the ratio $p_{2}/p_{1}$ must tend to $+\infty$ as $D_{21}\rightarrow +\infty$, since $D_{12}+\gamma \max\{p_{1}-p_{2}, 0\}$ is bounded independent of $D_{21}$ (Dividing the equilibria equation for $p_{1}$ by $D_{21}$ shows that $p_{1}\rightarrow 0$ as $D_{21}\rightarrow +\infty$ were there any such solutions).

However, if we add the equations for $p_{1}$ and $p_{2}$ we get

$0=(\mu s_{1}-d)p_{1}-\alpha p_{1}^2+(\mu s_{2}-d)p_{2}-\alpha p_{2}^2.$

Thus

$0=(\mu s_{1}-d)-\alpha p_{1}+(\mu s_{2}-d)\frac{p_{2}}{p_{1}}-\alpha\frac{p_{2}^2}{p_{1}}$

which implies

$[d-\mu s_{2}+\alpha p_{2}]\left(\frac{p_{2}}{p_{1}}\right)=\mu s_{1}-d-\alpha p_{1}.$

Recall that $\mu\beta_{1}-d>0$ and $\mu s_{2}-d\leq\mu\beta_{2}-d<0$. Consequently $d-\mu s_{2}+\alpha\beta_{2}>0$ and we get

$\frac{p_{2}}{p_{1}}=\frac{\mu s_{1}-d-\alpha p_{1}}{d-\mu s_{2}+\alpha p_{2}}\leq \frac{\mu\beta_{1}-d}{d-\mu\beta_{2}}$

independent of $D_{21}$. Consequently, (7) can have no equilibrium with $p_{1}$ and $p_{2}$ positive for sufficiently large $D_{21}$.

Propositions 1 and 2 and Theorem 4.1 extend to the analogues of (7) for an arbitrary number of micro habitat patches, so that permanence or uniform persistence is determined by the stable manifold of $\{(0, \beta_{1}, 0, \beta_{2},..., 0, \beta_{n})\}$. Indeed, we can readily obtain that the system is permanent so long as $\mu\beta_{i}-d>0$, $\mu\beta_{i}-d-D_{i-1, i}-D_{i+1,i}>0$ for some $i$ (with appropriate analogues if $i=1$ or $n$).

So the situation of primary interest is when $\mu\beta_{i}-d>0$ for one value of $i$ and negative for all others. For the sake of specificity, assume $\mu\beta_{2}-d>0$ and $\mu\beta_{i}-d<0$ for $i\neq 2$. Assuming a componentwise positive equilibrium, we may argue as in Proposition 3 that

$(D_{21}+\gamma \max\{p_{2}-p_{1}, 0\})\left(\frac{p_{1}}{p_{2}}\right)+(D_{23}+\gamma \max\{p_{2}-p_{3}, 0\})\left(\frac{p_{3}}{p_{2}}\right)\\ =D_{12}+D_{32}-(\mu s_{2}-d)+\alpha p_{2}.$

Consequently, if $D_{12}+D_{32}$ becomes large we get that $p_{1}/p_{2}$ or $p_{3}/p_{2}$ must become large. On the other hand, adding all the $p$ equations yields

$0=\sum\limits_{i=1}^{n}(\mu s_{i}-d)p_{i}-\alpha(\sum\limits_{i=1}^{n}p_{i}^2)$

which we can re-write as

$\sum\limits_{i\neq 2}(d-\mu s_{i}+\alpha p_{i})\left(\frac{p_{i}}{p_{2}}\right)=(\mu s_{2}-d-\alpha p_{2})$

Since $d-\mu s_{i}+\alpha p_{i}>0$ then for each $i\neq 2$, we get

$\frac{p_{i}}{p_{2}}\leq \frac{\mu s_{2}-d-\alpha p_{2}}{d-\mu s_{i}+\alpha p_{i}}\leq \frac{\mu\beta_{2}-d}{d-\mu\beta_{i}}.$

So there can not be permanence of the system if diffusion from patch 2 is too large. Again, the interpretation is that patch 2 effectively becomes a sink instead of a source if $D_{12}+D_{32}$ is too large.

---

5. Analysis Ⅱ: Equal diffusion rates

In this section we focus on the case where the diffusion rates of the bacteria are everywhere equal. Our purpose is to highlight further the parallel between the role of diffusion in our model as compared to its role in a reaction-diffusion analogue, such as [7]. Again, we consider the case where exactly one patch is a source; i.e., where $\mu\beta_{i}-d>0$ for exactly one $i$. In parallel to [7], we will also assume that $\sum_{i=1}^n (\mu\beta_{i}-d)<0$. The results from the previous section when $\gamma>0$ carry over when $D$ (the common diffusion rate) is small and may be framed as asserting the persistence of the bacteria in all MHP's in that case. Here we will take $\gamma=0$ and consider the predictions of the model as $D$ varies from $0$ to $\infty$. What we observe is that there is a critical value $\hat D$ of $D$ so that we get a prediction of persistence of bacteria in all MHP's when $0<D<\hat D$ and extinction when $D>\hat D$.

To this end, for the sake of specificity, we take $\mu\beta_{1}-d>0$ and $\mu\beta_{i}-d<0$ for $i=2,...,n$. Let $a_{i}=|\mu\beta_{i}-d|$ and assume that $a_{1}-a_{2}-a_{3}-...-a_{n}<0$. The predictions of the analogue to Theorem 4.1 can be discerned via the eigenvalues of the matrix

$A_1(D)= \begin{pmatrix} a_1-D & D & 0 & \ldots & 0 & 0 & 0\cr D & -a_2-2D & D & \ldots & 0 & 0 & 0\cr & \vdots & & \ddots & & \vdots & \cr 0 & 0 & 0 & \ldots & D & -a_{n-1}-2D & D \cr 0 & 0 & 0 & \ldots & 0 & D & -a_{n}-D \cr \end{pmatrix}$

Since $A_1(D)$ is a real symmetric matrix, all of its eigenvalues are real. The main result of this section is the following

Theorem 5.1. There is a unique positive number $\hat D$ such that for $0\le D<\hat D$ the matrix $A_1(D)$ has only one positive eigenvalue. The multiplicity of this eigenvalue is one and the other eigenvalues are all negative. For $D>\hat D$ the eigenvalues of $A_1(D)$ are all negative.

Remark 1. Since $\gamma=0$, the Hartman-Grobman Theorem implies that ${{W}^{s}}(\{(0,{{\beta }_{1}},\cdots ,0,{{\beta }_{n}})\})\cap {{X}_{1}}\ne \varnothing$ when $D>\hat D$.

Before proving this result, we introduce some notation and make some preliminary observations. We transform $A_1(D)$ by performing the following sequence of row operations: starting with row $i=1$ and ending with row $n-1$, add row $i$ to row $i+1$ which yields the matrix

$B_1(D)= \begin{pmatrix} a_1-D & D & 0 & \ldots & 0 & 0 & 0\cr a_1 & -a_2-D & D & \ldots & 0 & 0 & 0\cr a_1 & -a_2 & -a_3-D & \ldots & 0 & 0 & 0\cr & \vdots & & \ddots & & \vdots & \cr a_1 & -a_2 & -a_3 & \ldots & -a_{n-2} & -a_{n-1}-D & D \cr a_1 & -a_2 & -a_3 & \ldots & -a_{n-2} & -a_{n-1} & -a_{n} \cr \end{pmatrix}$

Next we perform the following sequence of column operations: starting with column $j=n$ and ending with column $2$, add column $j$ to column $j-1$. This yields the matrix

$C_1(D)=\begin{pmatrix} a_1 & D & 0 & \ldots & 0 & 0\cr a_1-a_2 & -a_2 & D & \ldots & 0 & 0\cr a_1-a_2-a_3 & -a_2-a_3 & -a_3 & \ldots & 0 & 0\cr & \vdots & & \ddots & \vdots & \cr c_{n-1,1} & c_{n-1,2} & c_{n-1,3} & \ldots & -a_{n-1} & D \cr c_{n,1} & c_{n,2} & c_{n,3} & \ldots & -a_{n-1}-a_{n} & -a_{n} \cr \end{pmatrix}$

where $c_{n,1}=a_1-a_2-a_3-\cdots-a_n$. Notice that the determinant is invariant under all these row and column operations. Thus we have the following

Lemma 5.2. The polynomial $p_1(D)=\det A_1(D)=\det C_1(D)$ is of degree $n-1$ in $D$ and

$p_1(D)=(-1)^{n-1}(a_1-a_2-a_3-\cdots-a_n)D^{n-1}+\cdots+(-1)^{n-1}a_1 a_2 \cdots a_n.$

Proof. The constant term of $p_1$ can easily be extracted by taking the determinant of the diagonal matrix $A_1(0)$. To extract the highest degree term, consider the determinant as the sum over all permutations of $\{1,\ldots,n\}$.

$\det {C_1}(D) = \sum\limits_{\sigma \in {S_n}} {{\mathop{\rm sgn}} } (\sigma )\prod\limits_{i = 1}^n {{c_{i,\sigma (i)}}}$

and notice that since $D$ does not appear in the last row, no term in this sum can be of degree higher than $n-1$. Further, the only term in this sum of degree $n-1$ is

$c_{1,2}c_{2,3}c_{3,4}\cdots c_{n-1,n}c_{n,1}=(a_1-a_2-a_3-\cdots-a_n)D^{n-1}$

and the corresponding permutation $\sigma$ is the $n$-cycle $(1,2,3,\ldots,n)$ which is even if $n$ is odd and odd if $n$ is even.

Notice that the lowest and highest degree terms of $p_1$ have opposite signs. Hence $p_1$ must have at least one root in $(0,\infty)$. Another thing to observe: by the Gershgorin Circle Theorem, the eigenvalues of $A_1(D)$ all lie in the union of intervals

$[a_1-2D,a_1]\cup\bigcup\limits_{j=2}^{n-1}[-a_j-4D,-a_j]\cup[-a_n-2D,-a_n].$

Notice that this implies that all positive eigenvalues must lie in the interval $(0,a_1]$. That leads to the following

Lemma 5.3. Let $D_1\le D_2$ be the smallest and largest positive roots of $p_1(D)$. Then for all $D$ in the interval $[0,D_1)$ the matrix $A_1(D)$ has exactly one positive eigenvalue. This eigenvalue is of multiplicity one and all other eigenvalues are negative. Further, for all $D$ in the interval $(D_2,\infty)$ all eigenvalues of $A_1(D)$ are all negative.

Remark 2. Theorem 5.1 follows once we show $D_1=D_2$.

Proof. The diagonal matrix $A_1(0)$ has exactly one positive eigenvalue with the rest negative. Since the eigenvalues of $A_1(D)$ are all real and depend continuously on $D$, as one increases $D$ from zero the positivity of one and the negativity of the rest of the eigenvalues can only change at a point where their product, $p_1(D)$, passes through zero. Thus $A_1(D)$ must have the desired properties for $0\le D<D_1$.

Similarly, for $D>D_2$ no eigenvalue of $A_1(D)$ can change sign. Since $\lim_{D\to\infty}p_1(D) = \infty$ if $n$ is even and $\lim_{D\to\infty}p_1(D)=-\infty$ if $n$ is odd, positive eigenvalues of $A_1(D)$ in $(D_2,\infty)$ must occur in pairs. Let $j$ be the number of pairs of these positive eigenvalues and let $\lambda_1(D),\ldots\lambda_{2j}(D)$ be the positive eigenvalues. Then by the Gershgorin Circle Theorem

$|p_1(D)|=|\lambda_1(D)\cdots\lambda_{2j}(D)\lambda_{2j+1}(D)\cdots\lambda_{n}(D)|\le a_1^{2j}\left(\max\limits_{2\le k \le n}(a_k)+4D\right)^{n-2j}.$

(12)

By Lemma 5.2, $p_1(D)$ is of degree $n-1$ in $D$ while the right-hand side of (12) is of degree $n-2j$ in $D$. But this can only be true for all $D$ sufficiently large if $j=0$.

As mentioned in the remark above, to complete the proof of Theorem 5.1, we must show that $p_1$ has no more than one positive root. This would follow if $p_1$ was monotone in $D$ on $(0,\infty)$, but there are counterexamples to this even when $n=2$. To this point in this section, a subscript on $A_1$ and $p_1$ may seem like an unnecessary bit of notation, but we now consider some other matrices.

For $1<k<n$ we define the principal minors of $A_1(D)$

$A_k(D)= \begin{pmatrix} -a_k-2D & D & 0 & \ldots & 0 & 0 & 0\cr D & -a_{k+1}-2D & D & \ldots & 0 & 0 & 0\cr & \vdots & & \ddots & & \vdots & \cr 0 & 0 & 0 & \ldots & D & -a_{n-1}-2D & D \cr 0 & 0 & 0 & \ldots & 0 & D & -a_{n}-D \cr \end{pmatrix}$

and take $A_n(D)=(-a_n-D)$. Further, we define the polynomials $p_k$ by

$p_k(D)=\det A_k(D), 1\le k\le n.$

As with $A_1(D)$, by the Gershgorin Circle Theorem, for $k>1$, the eigenvalues of the symmetric matrices $A_k(D)$ all lie in the union of intervals

$\bigcup\limits_{j = k}^{n - 1} {[ - {a_j} - 4D, - {a_j}] \cup [ - {a_n} - 2D, - {a_n}]}$

and thus the eigenvalues of $A_k(D)$ are all negative. Hence $A_k(D)$ is negative definite (i.e. $-A_k(D)$ is positive definite). Finally it is convenient to introduce symmetric positive definite matrices ${\hat A}_k(D)$, the negative of $A_k(D)$ with the diagonal reversed. That is, for $1<k<n$

${\hat A}_k(D)=\begin{pmatrix} a_n+D & -D & 0 & \ldots & 0 & 0 & 0\cr -D & a_{n-1}+2D & -D & \ldots & 0 & 0 & 0\cr & \vdots & & \ddots & & \vdots & \cr 0 & 0 & 0 & \ldots & -D & a_{k+1}+2D & -D \cr 0 & 0 & 0 & \ldots & 0 & -D & a_{k}+2D \cr \end{pmatrix}$

and ${\hat A}_n(D)=-A_n(D)=(a_n+D)$. Since flipping the diagonal of $A_k$ can be accomplished with an even total number of column and row swaps, we have

$\det(A_k(D))=(-1)^{n-k+1}\det({\hat A}_k(D))$

Since the matrices ${\hat A}_k$ are positive definite, they each have a Cholesky factorization $R_k$. That is, there exist upper triangular matrices $R_k$ with positive diagonal entries such that

$R_k^T\,R_k={\hat A}_k.$

Starting in the upper left corner, we let $r_n$ be the (only) entry of $R_n$ and hence

$r_n^2=a_n+D.$

(13)

Then for $2<k\le n$, by partitioning

$R_{k-1}=\begin{pmatrix}{\bar R}_{k} & v_{k-1} \cr 0 & r_{k-1}\cr\end{pmatrix}$

we see that

${\hat A}_{k-1} = R_{k-1}^T R_{k-1} = \begin{pmatrix}{\bar R}_{k}^T\,{\bar R}_{k} & {\bar R}_{k}^T\,v_{k-1}\cr v_{k-1}^T\,{\bar R}_{k} & v_{k-1}^T v_{k-1}+r_{k-1}^2\cr\end{pmatrix} = \begin{pmatrix}{\hat A}_{k} & b \cr b^T & a_{k-1}+2D \cr\end{pmatrix}$

where

$b^T = \left[\begin{matrix}0 & 0 & \ldots & 0 & -D\end{matrix}\right].$

Thus we see that in fact ${\bar R}_k=R_k$, that

$v_{k-1}^T = \left[\begin{matrix}0 & 0 & \ldots & 0 & -D/r_k\end{matrix}\right]$

and hence also that

$D^2/r_k^2+r_{k-1}^2=a_{k-1}+2D.$

As mentioned above, $p_1$ may not be strictly monotone for $D>0$, but by expanding $p_1$ along the first column of $A_1$ we see

$p_1(D)=(a_1-D)p_2(D)-D^2p_3(D).$

Since $\hat A_2(D)$ is positive definite for $D\ge0$, $p_2(D)$ is never zero there. We next consider the quotient

$Q(D)=\frac{p_1(D)}{p_2(D)} = a_1 - D - \frac{D^2\,p_3(D)}{p_2(D)}$

(14)

in our final

Lemma 5.4. The function $Q$ is strictly decreasing on $[0,\infty)$ with $Q(0)=a_1$and

$\lim\limits_{D\to\infty}Q(D)=a_1-a_2-a_3-\cdots-a_n.$

In particular, $Q$ has a unique root $\hat D$ in $[0,\infty)$.

Remark 3. Since the roots of $p_1$ and $Q$ coincide for $D>0$, this $\hat D$ is the same as in Theorem 5.1 and is also $D_1$ and $D_2$ in Lemma 5.3. Thus Theorem 5.1 follows immediately from Lemma 5.4.

Proof. From (14) we consider

$\begin{array}{l} {f_2}(D) = - \frac{{{D^2}{\mkern 1mu} {p_3}(D)}}{{{p_2}(D)}} - D\\ \quad \quad \quad = - \frac{{{D^2}{\mkern 1mu} \det ({A_3}(D))}}{{\det ({A_2}(D))}} - D\\ \quad \quad \quad = \frac{{{D^2}{\mkern 1mu} \det ({{\hat A}_3}(D))}}{{\det ({{\hat A}_2}(D))}} - D\\ \quad \quad \quad = \frac{{{D^2}{\mkern 1mu} \det {{({R_3})}^2}}}{{\det {{({R_2})}^2}}} - D\\ \quad \quad \quad = \frac{{{D^2}{\mkern 1mu} \det {{({R_3})}^2}}}{{\det {{({R_3})}^2}r_2^2}} - D\\ \quad \quad \quad = \frac{{{D^2}}}{{r_2^2}} - D. \end{array}$

Recall from (13) that $r_n^2=a_n+D$. We use induction starting with

$f_n(D) = \frac{D^2}{r_n^2}-D = \frac{D^2}{D+a_n}-D = \frac{-D\,a_n}{D+a_n}.$

Notice that

${f_{n'}}(D) = \frac{{ - a_n^2}}{{{{(D + {a_n})}^2}}} < 0,\quad \text{for}\;D \ge 0$

Further, $f_n(0)=0$ and $f_n(D)\to-a_n$ as $D\to\infty$.

Assume that for $2<k\le n$

$f_k(D)=\frac{D^2}{r_k^2}-D$

satisfies

${f_k}(D) < 0\;{\rm{for}}\;D > 0,$

(15)

${f_{k'}}(D) < 0\;{\rm{for}}\;D \ge 0,\;{\rm{and}}$

(16)

${f_k}(D) \to - {a_k} - {a_{k + 1}} - \cdots - {a_n}\;{\rm{as}}\;D \to \infty .$

(17)

Notice

$\begin{array}{l} {f_{k - 1}}(D) = \frac{{{D^2}}}{{r_{k - 1}^2}} - D\\ \quad \quad \quad \;\; = \frac{{{D^2}}}{{{a_{k - 1}} + 2D - {D^2}/r_k^2}} - D\\ \quad \quad \quad \;\; = \frac{{{D^2}}}{{{a_{k - 1}} + D - ({D^2}/r_k^2 - D)}} - D\\ \quad \quad \quad \;\; = \frac{{{D^2}}}{{{a_{k - 1}} + D - {f_k}(D)}} - D\\ \quad \quad \quad \;\; = \frac{{{D^2} - D({a_{k - 1}} + D - {f_k}(D))}}{{{a_{k - 1}} + D - {f_k}(D)}}\\ \quad \quad \quad \;\; = \frac{{ - D({a_{k - 1}} - {f_k}(D))}}{{D + ({a_{k - 1}} - {f_k}(D)).}} \end{array}$

Hence

$\begin{array}{l} {f_{k - {1^\prime }}}(D){\rm{ = }}\{ [ - ({a_{k - 1}} - {f_k}(D)) + D{f_{k'}}(D)][D + ({a_{k - 1}} - {f_k}(D))]{\rm{ = }}\\ \quad \quad \quad \quad \quad \quad + D({a_{k - 1}} - {f_k}(D))(1 - {f_{k'}}(D))\} /{(D + ({a_{k - 1}} - {f_k}(D)))^2}\\ \quad \quad \quad {\rm{ = }}\{ - D({a_{k - 1}} - {f_k}(D)) - {({a_{k - 1}} - {f_k}(D))^2} + {D^2}{f_{k'}}(D)\\ \quad \quad \quad \quad \quad \quad + D{f_{k'}}(D)({a_{k - 1}} - {f_k}(D)) + D({a_{k - 1}} - {f_k}(D))\\ \quad \quad \quad \quad \quad \quad - D{f_{k'}}(D)({a_{k - 1}} - {f_k}(D))\} /{(D + ({a_{k - 1}} - {f_k}(D)))^2}\\ \quad \quad \quad {\rm{ = }}\{ - {({a_{k - 1}} - {f_k}(D))^2} + {D^2}{f_{k'}}(D)\} /{(D + ({a_{k - 1}} - {f_k}(D)))^2} \end{array}$

which is negative for $D\ge0$. Since $f_{k-1}(0)=0$, $f_{k-1}(D)$ must be negative for $D>0$. Thus (15) and (16) of the induction hypothesis hold for $f_{k-1}$. Finally

$\mathop {\lim }\limits_{D \to \infty } {f_{k - 1}}(D) = \mathop {\lim }\limits_{D \to \infty } \frac{{ - {a_{k - 1}} + {f_k}(D)}}{{1 + ({a_{k - 1}} - {f_k}(D))/D.}} = - {a_{k - 1}} - {a_k} - \cdots - {a_n},$

which completes the induction step. Notice that

$Q(D)=a_1+f_2(D).$

and Lemma 5.4 follows. As mentioned above, this proves Theorem 5.1.

---

6. Numerical investigations

We have shown in the preceding sections that if there is a single source patch, the model predicts persistence of the bacteria in the micro-habitat patch array so long as the diffusion rates are not too large. The model is spatially implicit, but the underlying assumption that the array is linear (i.e. one must pass through patch $i$ to get to patch $(i+1)$ or patch $(i-1)$) enables us to think of the number of patch spaces between two patches as a surrogate for distance. In this section, we want to explore the asymptotic spatial arrangement of the population equilibrium distribution. It is certainly reasonable to expect that the population value in an MHP decreases as the bacteria is further away from a source patch, since it passes through MHP's that are net unfavorable. But it is interesting to consider how fast it decreases, and how is this impacted by the variation in net unfavorability among the sink patches and by the strength of the bacterial tendency to self-aggregate. In this section we explore numerically the input that various features of the model (6) and its analogues have upon the predictions of the models. We used the ode45 function in MATLAB to integrate numerically in the model. We begin in Experiment $1$ with the role of self-aggregation of bacteria, measured through the parameter $\gamma$. We illustrate the role of bacterial self-aggregation in (6) with 5 microhabitat patches (MHP's). In $(a)-(d)$ of Figure $1$, all parameters except $\gamma$ are fixed. Here $D_{ij}=D=0.01$, $\beta_{1}=0.2$, $\beta_{2}=0.2$, $\beta_{3}=0.8$, $\beta_{4}=0.02$, $\beta_{5}=0.38$, $\mu=0.15$, $\varepsilon=1.2$, $\lambda=0.004$, $d=0.06$ and $\alpha=0.002$. Only $\mu\beta_{3}-d>0$. $\gamma$ varies across $(a)-(d)$ with $\gamma=0$ $(a)$, $\gamma=4$ $(b)$, $\gamma=10$ $(c)$, $\gamma=50$ $(d)$. We see that increasing $\gamma$ leads to an equilibrium distribution where the abundance of the bacteria is increasingly concentrated in the source patch, patch 3. Consequently we can see that self-aggregation works to concentrate bacteria in the favorable habitat. See Figure 1.

Experiment 2 considers a case where there are 7 microhabitat patches. The middle patch, patch $4$, is favorable, meaning that $\mu\beta_{4}-d>0$, while $\mu\beta_{i}-d$ is negative for all other patches. Here $D_{ij}=D=0.01$, $\mu=0.15$, $\varepsilon=1.2$, $\lambda=0.004$, $d=0.06$, $\gamma=4$ and $\alpha$ is now $0.00002$. We explore the interplay between implicit distance from the unique source patch and relative disparity among maximal substrate input values. Here the $\beta_{i}$'s are ordered with $\beta_{1}<\beta_{2}<\beta_{3}<\beta_{4}>\beta_{5}>\beta_{6}>\beta_{7}$. We hold $\beta_{1}$, $\beta_{2}$, $\beta_{4}$, $\beta_{6}$ and $\beta_{7}$ fixed with $\beta_{1}=0.1$, $\beta_{2}=0.12$, $\beta_{4}=0.6$, $\beta_{6}=0.04$, $\beta_{7}=0.02$ so that $\beta_{1}=5\beta_{7}$ and $\beta_{2}=3\beta_{6}$. We initially have $\beta_{3}=0.3$ and $\beta_{5}=0.14$. In this instance we get $p_{1}>p_{7}$, $p_{2}>p_{6}$ and $p_{3}>p_{5}$, much as we would expect. We gradually reduce $\beta_{3}$ and increase $\beta_{5}$ till $\beta_{3}=0.14$ and $\beta_{5}=0.3$. When we have $\beta_{3}=0.2$ and $\beta_{5}=0.24$ so that $\beta_{3}$ is now less than $\beta_{5}$, we get $p_{3}<p_{5}$. However, at this point $p_{2}$ remains larger than $p_{6}$. By reducing $\beta_{3}$ to $0.15$ and increasing $\beta_{5}$ to $0.29$, we get that $p_{3}<p_{5}$ and also $p_{2}<p_{6}$, reflecting having sufficient disparity between $p_{3}$ and $p_{5}$. At this point, $p_{1}$ remains larger than $p_{7}$, reflecting that $\beta_{1}=5\beta_{7}$. Such remains the case when $\beta_{3}=0.14$ and $\beta_{5}=0.3$. If we now keep $\beta_{3}$ at $0.14$, but increase $\beta_{5}$ to $0.39$ (just below the breakeven point $0.4$ for $\mu\beta-d$), we get $p_{1}<p_{7}$. See Table 1.

---

7. Conclusions

This paper has been inspired by the nano-ecology experiments on the bacterium Escherichia coli by Keymer et al. described in [6] and the subsequent chemotactic reaction-diffusion model developed and analyzed in [3] which was prompted by the experiments. Our main aim was to model the system via an island-patch or discrete-diffusion system so as to capture the discrete nature of the micro-habitat patches (MHP's) within the overall array of patches and corridors. We also modified the formulation of logistic growth in the within-patch model so as to allow for net negative growth rates in individual patches.

The linear nature of the array of MHP's results in a tri-diagonal system in which the bacteria must pass through patch $i$ in order to transit from patch $i-1$ to patch $i+1$ or vice versa. Such a highly connected system leads to a discrete diffusion model which is irreducible (with or without the extenuating effect of chemotactic aggregation). The impact of irreducibility here is precisely analogous to that of the maximum principle in a diffusive model in a continuous space setting.

The prediction of the model is either that the bacteria persist in all MHP's or that they tend toward extinction in all patches. Moreover, one may use acyclicity results from persistence theory to see that which alternative obtains is determined by whether the stable set of the equilibrium with the bacteria absent contains any fully nontrivial initial configuration of the system (as in Theorem 4.1). Such is a consequence of the Acyclicity Theorem of persistence theory via the results of [9].

The model exhibits strong source-sink dynamics when resource flow into some MHP's is set low enough so that a bacteria population is not sustainable in such patches in isolation. In such instances, diffusive dispersal may serve as a rescuing mechanism. Indeed, if there is a single patch in which the bacteria can survive in isolation (what we term a favorable patch), a slow rate of diffusion from the favorable patch leads to coexistence in all patches. However, if the rate of diffusion from the favorable patch is too high relative to diffusion into it from adjacent patches, the rescue effect is insufficient and the bacteria tend to extinction in the system. Such is the case whether or not there is bacterial self-aggregation in the system.

Of course, such a disparity in dispersal is not possible in the special but natural case when the diffusion rates are the same in all patches. Here we consider the case when patch 1 is the sole favorable patch and the overall habitat is unfavorable in the sense that the sum over all patches of net growth rates is negative. This assumption is analogous to the assumption that the integral of the growth rate is negative in [7]. In this case, when there is no bacterial self-aggregation, persistence is equivalent to the instability of the bacteria absent equilibrium. We show in Theorem 2 that there is a unique positive threshold value of the diffusion rate D so that the bacteria absent equilibrium is unstable for diffusion rates below the threshold and stable for values above the threshold. Such is the case even though the determinant of the relevant Jacobi matrix is not monotonic as a function of the diffusion rate.

Based on our numerical experiments, the effect of bacterial self-aggregation appears to be to concentrate the population in favorable MHP's. In our two experiments, we consider the situation where we have 5 and 7 patches wherein the middle patch is favorable while the overall environment is net unfavorable in the sense we have described. The long term population density is positive in all patches but trails off when one moves away from the favorable patch. As the self-aggregation parameter is increased in Experiment $1$ the disparity between the density in the favorable patch and the density elsewhere becomes more and more pronounced. Experiment 2 suggests that the interplay of patch "distance" and maximal substrate rates is nuanced. In the experiment, we kept maximal substrate rates constant in patches $1$, $2$, $4$, $6$ and $7$ and varied them in patches $3$ and $5$, with $\beta_{1}<\beta_{2}<\beta_{3}<\beta_{4}>\beta_{5}>\beta_{6}>\beta_{7}$, with $\beta_{1}=5\beta_{7}$ and $\beta_{2}=3\beta_{6}$, and with only $\mu\beta_{4}-d>0$. When $\beta_{3}$ is considerably larger than $\beta_{5}$, we found the equilibrium populations (the $p_{i}$'s) match the relations among the maximal substrate rates (the $\beta_{i}$'s). As $\beta_{3}$ decreases and $\beta_{5}$ increases, we first reverse the equilibrium sizes of $p_{3}$ and $p_{5}$. If the disparity between $\beta_{3}$ and $\beta_{5}$ is not too large, it remains the case that $p_{2}>p_{6}$ and $p_{1}>p_{7}$, suggesting that the disparity in maximal substrate rates ($\beta_{1}=5\beta_{7}$ and $\beta_{2}=3\beta_{6}$) is outweighing the feed from the inner most unfavorable patches ($p_{3}$ and $p_{5}$). Once $\beta_{3}$ is small enough relative to $\beta_{5}$, we find that $p_{2}<p_{6}$ even though $\beta_{2}=3\beta_{6}$. If $\beta_{3}$ is still further smaller than $\beta_{5}$, $p_{1}<p_{7}$, even though $\beta_{1}=5\beta_{7}$.

---