A tridiagonal patch model of bacteria inhabiting a nanofabricated landscape

In this paper we employ a discrete-diffusion modeling framework to examine a system inspired by the nano-ecology experiments on the bacterium Escherichia coli reported upon in Keymer et al. (2006). In these experiments, the bacteria inhabit a linear array of 85 ``microhabitat patches (MHP's)", linked by comparatively thinner corridors through which bacteria may pass between adjacent MHP's. Each MHP is connected to its own source of nutrient substrate, which flows into the MHP at a rate that can be controlled in the experiment. Logistic dynamics are assumed within each MHP, and nutrient substrate flow determines the prediction of the within MHP dynamics in the absence of bacteria dispersal between patches. Patches where the substrate flow rate is sufficiently high sustain the bacteria in the absence of between patch movement and may be regarded as sources, while those with insufficient substrate flow lead to the extinction of the bacteria in the within patch environment and may be regarded as sinks. We examine the role of dispersal in determining the predictions of the model under source-sink dynamics.


(Communicated by Yuan Lou)
Abstract. In this paper we employ a discrete-diffusion modeling framework to examine a system inspired by the nano-ecology experiments on the bacterium Escherichia coli reported upon in Keymer et al. (2006). In these experiments, the bacteria inhabit a linear array of 85 "microhabitat patches (MHP's)", linked by comparatively thinner corridors through which bacteria may pass between adjacent MHP's. Each MHP is connected to its own source of nutrient substrate, which flows into the MHP at a rate that can be controlled in the experiment. Logistic dynamics are assumed within each MHP, and nutrient substrate flow determines the prediction of the within MHP dynamics in the absence of bacteria dispersal between patches. Patches where the substrate flow rate is sufficiently high sustain the bacteria in the absence of between patch movement and may be regarded as sources, while those with insufficient substrate flow lead to the extinction of the bacteria in the within patch environment and may be regarded as sinks. We examine the role of dispersal in determining the predictions of the model under source-sink dynamics.
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µm (length) ×100µm (width) ×30µm (depth). By comparison, a corridor connecting two adjacent MHP's is 50µm (length) ×5µm (width) ×30µm (depth). Each MHP is separately attached via two feeding channels to the nutrient source. Each channel has 5 "nanoslits" of dimensions 20µm (length) ×15µm (width) ×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 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 where α > 0. When r > 0, (2) is equivalent to (1) with K = r/α. When r ≤ 0, x(t) tends to zero as t → ∞ 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 bacterianutrient 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/α 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 where p is the density of the bacteria, s is the density of the substrate that the bacteria feed upon, µ is the bacterial growth rate under the maximal sustainable substrate level and d is bacterial mortality. K represents the carrying capacity of the habitat, λ is the flow rate through the feeding corridor, and ε is a conversion factor from the consumption of substrate to the birth of bacteria. In this formulation, ε > 1. As discussed in the introduction, we modify (3) to (4) below to allow for the possibility that µ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 β ∈ (0, 1] for this purpose. Consequently, our single patch model is expressed as 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 which converges over time to β. If µβ − d < 0, then the bacterial density is a lower solution to for t 0 and consequently (p, s) → (0, β) as t → ∞. A modified argument holds when µβ − d = 0. So now suppose µβ − d > 0. An equilibrium (p * , s * ) with p * = 0 satisfies so that s * must satisfy so that is the only positive root of (5). Setting f (s) = εµ 2 Ks 2 + (λ − εµKd)s − λβ, we have f (0) = −λβ while f (β) = εµKβ(µβ − d) > 0 so that s * < β. Additionally, f (d/µ) = λ(d/µ − β) < 0 so that d/µ < s * < β and 0 < p * < K(µβ − 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 In particular, (p * , s * ) lies within the rectangle. At (p * , s * ) the Jacobian matrix simplifies to 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 in the interior of D. Consequently, we have by Dulac's Criterion and Poincaré-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 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 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 where i = 2, . . . , 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), γ denotes the strength of chemotactic bacterial self-aggregation, and ν the strength of bacterial cueing upon relative substrate concentrations. The single patch parameters µ,d,λ and ε have the same meanings as in the preceding section, as does also β i , except β 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 α in the bacterial logistic dynamics. The parameters γ, ν, µ, d, ε and α reflect bacterial traits which can be reasonably considered as independent of patch. Similar reasoning applies to λ relative to the substrate.
In the analysis that follows, we will set ν = 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 is bounded by νβ i+1 p i . Hence it is relatively insignificant as a dispersal term compared to the linear movement terms if ν 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 ≡ 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 I: General diffusion rates. We will give persistence results for the n-patch analogue of (6) with ν = 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 dp Here we are interested in the solution flow for (7) on the set X = {(p 1 , s 1 , p 2 , s 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 so we start our discussion with a verification of this fact.
is a solution of (7) with p i (0) > 0 and 0 < s i (0) < β 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 In such case, either s i (t 1 ) = β i or 0 for at least one i or p i (t 1 ) = 0 for at least one i.

So consider
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 (7) is not forward invariant on ∂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 d dt Consequently, given any σ > 0, there is a T = T (p 1 (0), p 2 (0)) such that for t ≥ 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 ν = 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 ≡ β i or p i ≡ 0 for i = 1 or i = 2. It is immediate that s i ≡ β i for i = 1 or 2 implies that p i ≡ 0. The tridiagonal structure of the model at the linear level forces p j ≡ 0 where j = i. Having p j ≡ 0 in turn implies that s j → β j as t → ∞. Consequently, the only element of w(∂X 1 ) is {(0, β 1 , 0, β 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, β 1 , 0, β 2 )}, denoted W s ({(0, β 1 , 0, β 2 )}) ∩ X 1 = ∅. So we have: Theorem 4.1. The system (7) is permanent (uniformly persistent) if and only if We now examine when Theorem 4.1 holds. We begin with some simple observations. Note from (7) that If µβ 1 − d < 0 and µβ 2 − d < 0, we have by (8) that p 1 + p 2 is a subsolution of the initial value problem Consequently, for any initial data with p i (0) ≥ 0, we get that p 1 → 0 and p 2 → 0 as t → ∞, which in turn implies that W s ({(0, β 1 , 0, β 2 )}) ∩ X 1 = ∅. On the other hand, suppose that µβ 1 − d > 0 and µβ 2 − d > 0 and that (p 1 (t), s 1 (t), p 2 (t), s 2 (t)) → (0, β 1 , 0, β 2 ) for some initial data (p 1 (0), s 1 (0), p 2 (0), s 2 (0)) with p i (0) ≥ 0 and p 1 (0) + p 2 (0) > 0. Then p i (t) > 0 and 0 < s i (t) < β i for all t > 0 and for t sufficiently large, say t ≥ T 1 , p 1 + p 2 is a supersolution of the initial value problem where c is positive and less than min{µβ 1 − d, µβ 2 − d}.
Consequently, given any σ ∈ (0, c/α), there is a T 2 > T 1 so that p 1 (t) + p 2 (t) > c/α − σ for all t ≥ T 2 . In particular, it is not possible for a solution of (7) starting in X 1 to converge to 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 µβ i − d is positive in one patch (meaning a favorable or source habitat) but negative in the other (meaning an unfavorable or sink habitat).
The Jacobi matrix at (0, β 1 , 0, β 2 ), J(0, β 1 , 0, β 2 ) is given by We will use σ to denote the eigenvalues of J(0, β 1 , 0, β 2 ) in (9). It is not difficult to see that σ = −λ is a double eigenvalue (reflecting the restorative tendency of the substrate when bacterial abundance is low) and the other eigenvalues of J(0, β 1 , 0, β 2 ) are the eigenvalues of the matrix A given by So the other two eigenvalues are given by 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 , ∆ < 0 while T r could be of either sign. When D 12 and D 21 are positive and small, ∆ remains negative. Consequently, J(0, β 1 , 0, β 2 ) admits one positive eigenvalue and we have that W s ({(0, β 1 , 0, β 2 )}) ∩ X 1 = ∅, 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 T r < 0 when D 12 = 0 = D 21 , it remains negative for all positive values of D 12 and D 21 . If T r > 0 when D 12 = 0 = D 21 , it will eventually become negative as D 21 increases. By re-writing ∆ and T r as it is easy to see that ∆ remains negative as D 21 increases until after the point at which T r becomes negative as D 21 increases. It is also easy to see from the above that ∆ 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, β 1 , 0, β 2 )})∩ X 1 = ∅. 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 γ > 0, so that we take bacterial selfaggregation into account. Theorem 4.1 remains valid. Suppose that D 21 is small enough so that µβ 1 − d − D 21 > 0. Let r 1 and r 2 be positive numbers so that and suppose there is a solution of (7) with p i (0) > 0, s i (0) ∈ (0, β i ), i = 1, 2 so that (p 1 (t), s 1 (t), p 2 (t), s 2 (t)) → (0, β 1 , 0, β 2 ) as t → ∞. Choose T > 0 so that for t ≥ T , Then Hence we obtain that for t ≥ T , dp 1 dt as t → ∞, we get that p 1 (t) > w for t T , a contradiction. So there can be no such orbit and thus So (7) remains permanent when γ > 0, D 12 > 0 and D 21 is positive and small. When γ = 0 and D 12 > 0 is fixed, (7) loses permanence when D 21 becomes large enough so that the determinant ∆ of A in (10) becomes positive. For γ > 0 and D 12 > 0 and fixed, we establish that (7) is not permanent for D 21 sufficiently large via the following result. Proof. Suppose (7) admits an equilibrium with p i > 0 and s i ∈ (0, β i ) for i = 1, 2. Then we have that For the time dependent problem, we have d dt so that p 1 and p 2 are ultimately bounded above by any number greater than [2(µβ 1 − d)/α], 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 +∞ as D 21 → +∞. 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 +∞ as D 21 → +∞, since D 12 + γ 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 → 0 as D 21 → +∞ were there any such solutions).
However, if we add the equations for p 1 and p 2 we get Recall that µβ 1 − d > 0 and µs 2 − d ≤ µβ 2 − d < 0. Consequently d − µs 2 + αβ 2 > 0 and we get 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, β 1 , 0, β 2 , ..., 0, β n )}. Indeed, we can readily obtain that the system is permanent so long as So the situation of primary interest is when µβ i − d > 0 for one value of i and negative for all others. For the sake of specificity, assume µβ 2 − d > 0 and µβ i − d < 0 for i = 2. Assuming a componentwise positive equilibrium, we may argue as in Proposition 3 that 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 Since d − µs i + αp i > 0 then for each i = 2, we get 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 II: 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 µβ i − d > 0 for exactly one i. In parallel to [7], we will also assume that The results from the previous section when γ > 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 γ = 0 and consider the predictions of the model as D varies from 0 to ∞. What we observe is that there is a critical valueD of D so that we get a prediction of persistence of bacteria in all MHP's when 0 < D <D and extinction when D >D.
To this end, for the sake of specificity, we take µβ 1 − d > 0 and µβ i − d < 0 for i = 2, ..., n. Let a i = |µβ 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 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 . . . −a n−2 −a n−1 − D D a 1 −a 2 −a 3 . . . −a n−2 −a n−1 −a n          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 n−1,1 c n−1,2 c n−1,3 . . . −a n−1 D c n,1 c n,2 c n,3 . . . −a n−1 − a n −a n where c n,1 = a 1 − a 2 − a 3 − · · · − a n . Notice that the determinant is invariant under all these row and column operations. Thus we have the following a n )D n−1 + · · · + (−1) n−1 a 1 a 2 · · · 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, . . . , n}.
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 · · · c n−1,n c n,1 = (a 1 − a 2 − a 3 − · · · − a n )D n−1 and the corresponding permutation σ is the n-cycle (1, 2, 3, . . . , 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, ∞). Another thing to observe: by the Gershgorin Circle Theorem, the eigenvalues of A 1 (D) all lie in the union of intervals Notice that this implies that all positive eigenvalues must lie in the interval (0, a 1 ]. That leads to the following (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, ∞), 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) and take A n (D) = (−a n − D). Further, we define the polynomials p k by That is, for 1 < k < n andÂ 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(Â k (D)) Since the matricesÂ 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 T k R k =Â k . Starting in the upper left corner, we let r n be the (only) entry of R n and hence r 2 n = a n + D. (13) Then for 2 < k ≤ n, by partitioning we see that Thus we see that in factR k = R k , that v T k−1 = 0 0 . . . 0 −D/r k and hence also that D 2 /r 2 k + r 2 k−1 = 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 SinceÂ 2 (D) is positive definite for D ≥ 0, p 2 (D) is never zero there. We next consider the quotient in our final Lemma 5.4. The function Q is strictly decreasing on [0, ∞) with Q(0) = a 1 and lim D→∞ Q(D) = a 1 − a 2 − a 3 − · · · − a n .

Remark 3.
Since the roots of p 1 and Q coincide for D > 0, thisD 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 Recall from (13) that r 2 n = a n + D. We use induction starting with Notice that f n (D) = −a 2 n (D + a n ) 2 < 0, for D ≥ 0. Further, f n (0) = 0 and f n (D) → −a n as D → ∞.
Assume that for 2 < k ≤ n Notice ).
Hence 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 β 1 < β 2 < β 3 < β 4 > β 5 > β 6 > β 7 , with β 1 = 5β 7 and β 2 = 3β 6 , and with only µβ 4 − d > 0. When β 3 is considerably larger than β 5 , we found the equilibrium populations (the p i 's) match the relations among the maximal substrate rates (the β i 's). As β 3 decreases and β 5 increases, we first reverse the equilibrium sizes of p 3 and p 5 . If the disparity between β 3 and β 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 (β 1 = 5β 7 and β 2 = 3β 6 ) is outweighing the feed from the inner most unfavorable patches (p 3 and p 5 ). Once β 3 is small enough relative to β 5 , we find that p 2 < p 6 even though β 2 = 3β 6 . If β 3 is still further smaller than β 5 , p 1 < p 7 , even though β 1 = 5β 7 .