Mathematical Analysis of Two Competing Cancer Cell Migration Mechanisms Driven by Interstitial Fluid Flow

Recent experimental work has revealed that interstitial fluid flow can mobilize two types of tumor cell migration mechanisms. One is a chemotactic-driven mechanism where chemokine (chemical component) bounded to the extracellular matrix (ECM) is released and skewed in the flow direction. This leads to higher chemical concentrations downstream which the tumor cells can sense and migrate toward. The other is a mechanism where the flowing fluid imposes a stress on the tumor cells which triggers them to go in the upstream direction. Researchers have suggested that these two migration modes possibly can play a role in metastatic behavior, i.e., the process where tumor cells are able to break loose from the primary tumor and move to nearby lymphatic vessels. In Waldeland and Evje (J Biomech 81:22–35, 2018), a mathematical cell–fluid model was put forward based on a mixture theory formulation. It was demonstrated that the model was able to capture the main characteristics of the two competing migration mechanisms. The objective of the current work is to seek deeper insight into certain qualitative aspects of these competing mechanisms by means of mathematical methods. For that purpose, we propose a simpler version of the cell–fluid model mentioned above but such that the two competing migration mechanisms are retained. An initial cell distribution in a one-dimensional slab is exposed to a constant fluid flow from one end to the other, consistent with the experimental setup. Then, we explore by means of analytical estimates the long-time behavior of the two competing migration mechanisms for two different scenarios: (i) when the initial cell volume fraction is low and (ii) when the initial cell volume fraction is high. In particular, it is demonstrated in a strict mathematical sense that for a sufficiently low initial cell volume fraction, the downstream migration dominates in the sense that the solution converges to a downstream-dominated steady state as time elapses. On the other hand, with a sufficiently high initial cell volume fraction, the upstream migration mechanism is the stronger in the sense that the solution converges to an upstream-dominated steady state. Communicated by Alain Goriely. Extended author information available on the last page of the article


Aggressive Cancer Cells and Fluid Flow
How and why is it so that aggressive cancer cells are able to detach from the primary tumor and migrate to nearby lymphatic vessels through which they can escape and give rise to formation of tumors at other places in the human body? This phenomenon of lymph node metastasis, which is a main reason why cancer becomes a deadly disease, has been recognized for a long time. However, the underlying mechanism by which malignant tumor cells leave the primary tumor site, invade the lymphatics and metastasize to lymph nodes is unclear (Shields et al. 2007;Polacheck et al. 2011). Many malignant tumors are associated with an elevated interstitial fluid pressure (IFP) caused by leaky blood vessels situated at the inside of the punctum. Lymphatic vessels normally adsorb this fluid and keep the IFP at a normal level. However, lymphatic vessels are often defective in the intratumoral region. This implies that the additional fluid oozes to the region outside the tumor periphery where it is adsorbed by collecting lymphatic vessels. It has been proposed that this elevated IF flow can be exploited by the tumor cells and has led researchers to systematically explore how tumor cells are sensitive to IF flow. In Shields et al. (2007), it was suggested that interstitial flow caused by lymphatic drainage directs tumor cell migration through chemotaxis. More precisely, the tumor cells utilize interstitial flow to create and amplify gradients in chemokine (a protein) and thus chemotact toward the adsorbing lymphatic vessels in a process termed autologous chemotaxis. Polacheck et al. Polacheck et al. (2011) extended the study by Shields et al. (2007), demonstrating that the IF velocity as well as the cell seeding density affected the migration direction. Experiments were conducted at two different seeding densities and at two different flow velocities. In particular, it was observed that for the low cell seeding density, culture tumor cells tended to migrate with the flow in accordance with the behavior reported in Shields et al. (2007). However, for the high cell seeding density, the migration was dominated by upstream migration.

A General Cell-Fluid-ECM Model
A rather general cell-fluid-ECM model was proposed in Waldeland and Evje (2018a) and further developed in Waldeland and Evje (2018b) and Evje and Waldeland (2019) to shed light on the above-mentioned competing cell migration mechanisms governed by interstitial fluid flow. A gently simplified version of the model, where we ignore certain details of the biochemical part by assuming that chemokine C is directly produced by the tumor cells instead of being released from ECM, takes the following form: The model, which bears similarity to the model studied in Evje and Wen (2018), accounts for the volume fraction α w of interstitial fluid (IF) and the volume fraction α c of cancer cells such that α c + α w = 1. In other words, the pore space is occupied by cancer cells and fluid and described my the two mass balance equations (1.1) 1,2 . The two different phases move with their own interstitial velocity, respectively, u w and u c . These are involved in the two momentum balance equations (1.1) 3,4 . The momentum balance for the IF given by (1.1) 4 reflects that the interstitial fluid pressure gradient ∇ P w is balanced by two interaction forces whose coefficients areζ w andζ . The first one reflects the resistance force felt by the fluid as it flows through the porous tissue, whereas the second reflects a drag force effect between the fluid and cells. Similarly, the momentum balance (1.1) 3 reflects that the cell phase pressure P w + P(α c ) + (C) differs from the IFP P w by two stress effects: P(α c ) is an increasing function which accounts for the effect that cells tend to move away from each other toward a region of lower cell volume fraction when they are densely packed to reduce the total cell phase pressure. (C) is a decreasing function which accounts for the cell's ability to create directed motion toward higher concentration of C (i.e., toward positive gradients in C) to reduce the overall pressure. Similarly,ζ c represents cell-ECM interaction andζ the cell-fluid drag. The last equation (1.1) 5 reflects that the chemokine concentration C is advected according to the fluid velocity field u w , in addition to a diffusive spreading, combined with production and consumption as described by the source term R C .
From the two momentum equations (1.1) 3,4 , we can compute explicit expressions for the cell and fluid velocity, respectively, u c and u w (Waldeland and Evje 2018a). The following expressions are found: Waldeland and Evje (2018b) and where the total velocity U T = α c u c + α w u w is determined from the equation U T = −λ T ∇ P w −λ c ∇( P + ) and the fact that ∇ ·U T = 0. We refer to Waldeland and Evje (2018a) for details. Model (1.1) then takes the more compact form with unknowns (α c , C): where α c u c is given by (1.2) 1 . Note thatf c (α c ) andĥ(α c ) given by (1.3) are direct functions of the specified fluid-ECM interactionζ w , cell-ECM interactionζ c and cell-fluid interactionζ . These correlations reflect essential information in what way tumor cells respond and relate to their microenvironment. Moreover, we note that there are three different mechanisms involved in (1.2) 1 : (i) the term U Tfc (α c ) represents a cell migration effect due to fluid stress; (ii)ĥ c (α c )∇( P(α c )) represents a diffusive cell-cell migration effect; and (iii)ĥ c (α c )∇( (C)) represents a chemotaxis migration effect. For typical correlations used forζ w ,ζ c andζ , the shape off c (α c ) andĥ(α c ) will be as shown in Fig. 1. The resulting cell migration behavior is shown in Figs. 2 and 3, respectively, for the case with an initial low cell volume fraction and a high initial cell volume fraction. The numerical examples illustrate the competition between downstream and upstream migration as a function of cell volume fraction.

A Toy Model with Competing Downstream and Upstream Migration
In order to obtain a model that is more amenable for mathematical investigations, without losing the key characteristics of the cell-fluid model (1.4), we make the following assumptions: (i) We use the approximation U T ≈ u w ≈ const.

Fig. 2
Competing tumor cell migration mechanisms for a cell aggregate with low volume fraction α c ≈ 0.1. a cell volume fraction α c . The downstream chemotactic-driven mechanism dominates. b The interstitial fluid pressure P w gradient. High pressure at x = 0 and low pressure at x = 1 give rise to fluid flow from left to right. c The three different cell migration components. d Chemokine (chemical component) whose concentration is skewed in the downstream direction The IF velocity u w typically is a 100-fold higher than the cell migration velocity u c (Polacheck et al. 2011;Waldeland and Evje 2018a, b), which in turn is largely dictated by the linear pressure curve seen in Figs. 2 and 3 (panel B). This justifies assumption (i). The constant diffusion coefficient in (ii) is standard. The choice off c (α c ) in (iii) accounts for the negative dip that gives rise to upstream migration for higher cell volume fraction α c , see Fig. 1 (left). The choice ofĥ(α c ) in (iv) is also consistent with the functional form ofĥ(α c ) in (1.3) 3 which amounts to a bell-shaped function starting and ending at 0, see Fig. 1 (right), combined with the fact that (C) is a decreasing function Waldeland and Evje (2018a, b). Finally, the choice of parameters and terms in (v) is standard.
With these assumptions and by replacing α c and C by u and v, respectively, we obtain the following simplified version of (1.4): (1.5)

Fig. 3
Competing tumor cell migration mechanisms. The situation is the same as in Fig. 2. The only difference is that we use the higher volume fraction α c ≈ 0.5. The upstream mechanism now dominates.
The reason for that can be seen from an inspection off c (α c ) in Fig. 1. A larger cell volume fraction α c means that a larger part of the downward dip is activated and therefore increases the impact from the upstream migration where f (u) and h(u) are given by where κ > 1 and λ > 1 are fixed parameters.

Analysis of Related Models for Chemotaxis-Fluid Interplay in the Literature
Understanding the interaction of chemotaxis systems with liquid environments has been the objective of a remarkably quickly growing literature during the past few years. Most analytical studies in this direction, however, focus on models addressing situations in which besides fluid-induced transport mechanisms, also certain buoyancydriven gravitational forcing of the considered fluid flows is relevant; especially due the fact that then the fluid velocity forms a genuine unknown in the model, such additional feedback effects evidently go along with a noticeably higher complexity in comparison with (1.5); therefore, the analysis of accordingly obtained chemotaxis-(Navier-)Stokes systems (Tuval et al. 2005) has yet been predominantly concerned with rather basic issues such as questions from existence and regularity theory (Duan et al. 2010;Winkler 2012Winkler , 2016Chae et al. 2014;Cao and Lankeit 2016;Kozono et al. 2016), and only few studies seem go beyond this by examining qualitative aspects such as large-time stabilization toward homogeneous equilibria (Lankeit 2016;Winkler 2014Winkler , 2017Winkler , 2019. Only for some more specialized and simplified variants involving suitably designed given fluid flows, more subtle findings on possible effects of fluid interaction, e.g., on taxis driven blowup or also certain spreading properties, are available (He and Tadmor 2019;Kiselev and Xu 2016;Kiselev and Ryzhik 2012).

Main Results I: Dominance of Downstream Migration in Sparsely Distributed Populations
In the first part of our analysis, we shall consider the fully no-flux-type initial-boundary value problem (1.7) in the interval = (0, L) with L > 0, with suitably small and appropriately regular initial data in the sense that ε > 0 is suitably small and that (1.8) Within this framework, the first of our main results states that indeed for appropriately small values of ε, this problem is classically solvable by functions which exhibit a certain tendency toward migration in the direction of the fluid flow, that is, toward large positive values of x, in the following sense.

Main Results II: Prevalence of Upstream Migration in Densely Populated Groups
In the second part of our study, we consider the parabolic system from (1.7) along with slightly different boundary conditions for the first solution component, which namely enforce the latter to attain the boundary value 1 on ∂ throughout evolution. Specifically, for ε > 0 we shall be concerned with the problem (1.11) under the assumptions that (1.12) Under these assumptions, we shall see that in sharp contrast to the above, the profiles of the deviations 1−u ε from the level 1 will, throughout arbitrarily large time intervals, to a considerable extent remain near functions that are more concentrated near x = 0 than near x = L.

Classical Solutions to (1.7) in Ä × (0, T) for Small "
In order to construct solutions to (1.7) by means of a convenient approximation involving homogeneous Neumann boundary conditions in the first solution component, following precedent works pursuing a similar idea, we fix a family ( x ∈ , (2.1) admit local classical solutions enjoying a handy extensibility criterion: Lemma 2.1 Let ε > 0 and j ∈ N. Then, there exist T ε j and a unique pair of nonnegative functions Proof All statements can be verified by straightforward adaptation of well-known arguments to the present context, either following precedents concerned with taxistype problems, such as e.g., Horstmann and Winkler (2005), or also directly resorting to general theory for abstract parabolic evolution problems (Amann 1989).
A first significant regularity information about these solutions, becoming important in our derivation of uniform bounds on u ε j in Lemma 2.3, can be obtained by conveniently transforming the second equation in (2.1) and then performing an essentially straightforward testing procedure. As . Now in view of the Hölder inequality, for verifying the statement of the lemma, it is sufficient to establish (2.4) for any fixed integer q ≥ 2, and in order to achieve this, we again use the Neumann-type structure of the boundary condition in (2.5) to see that for any such q, Here, observing that by (2.6), (2.7), Young's inequality and the fact that c 1 (T ) ≥ 1, we can estimate we see that once more due to Young's inequality, , which in conjunction with (2.7) readily entails (2.4).
Thus, in particular, having at hand some information on integrability of the taxic gradient in (2.1), by making essential use of the presence of homogeneous Neumann boundary conditions for u ε j , we can invoke smoothing estimates for the Neumann heat semigroup to assert a favorable uniform a priori bound for the first solution component, up to an arbitrary fixed time.
Given T > 0, on employing Lemma 2.2, we can find c 1 (T ) > 0 such that for all ε > 0 and j ∈ N, (2.9) where again T ε j := min{T , T ε j }. We furthermore recall a well-known smoothing property of the Neumann heat semigroup (e τ ) τ ≥0 on (Fujie et al. 2016 (2.10) As κ > 1, it thereafter becomes possible to firstly pick δ(T ) > 0 small enough such that δ(T ) ≤ 1 and and then choose ε down (T ) > 0 in such a way that ( 2.13) We now fix ε ∈ (0, ε down (T )), and we claim that then for each j ∈ N, In fact, if this was false, then again by continuity of u ε j , we would have (2.14) To see that this is impossible, we represent u ε j according to a Duhamel formula associated with the first equation in (2.1) and apply the maximum principle as well as (2.10) to infer that for all t ∈ (0, T ε j ], and since by the Cauchy-Schwarz inequality, (2.3), (2.9) and the inequality δ(T ) ≤ 1, we moreover have this entails that thanks to (2.11), (2.12) and (2.13), When evaluated at t = T ε j , this contradicts (2.14) and thereby shows that indeed T ε j = T ε j . As thus u ε j (·, t) L ∞ ( ) ≤ δ(T ) ≤ 1 for all t ∈ (0, T ε j ), in view of (2.9) and (2.2) this furthermore implies that we must have T ε j > T ε j and that (2.8) holds.
On the basis of the latter, straightforward application of parabolic regularity theory, followed by suitable compactness arguments, enables us to construct a solution of (1.7) in × (0, T ) as a limit of solutions to (2.1), provided that ε < ε down (T ).

2.16)
This solution can be obtained as limits of the solutions to (2.1) in the sense that there exists a sequence ( j k ) k∈N ⊂ N such that as k → ∞, we have j k → ∞, u ε j k → u ε Proof Relying on Lemma 2.3 and a series of straightforward parabolic bootstrap arguments, thanks to the assumed limit behavior of (ζ j ) j∈N as j → ∞, the part concerning existence and approximation can be seen by following a type of reasoning well-established in contexts of taxis problems involving no-flux boundary conditions different from homogeneous Neumann data; as concise derivations can be found in quite an elaborate manner in the literature on closely related problems, we may refrain from giving details here, and rather refer to, for example, Cao and Lankeit (2016) (see also Li et al. (2015) for a precedent). The properties (2.15) and (2.16) can thereupon easily be obtained on taking j → ∞ in (2.3) and (2.8).
With regard to the rescaled version of u ε addressed in the finally intended estimate (1.10), this result can be rephrased as follows.

Existence, Uniqueness and Stabilization in a Formally Obtained Limit Problem
Motivated by formally taking ε 0 in the reformulation (2.17) of (1.7), in this section we shall analyze the behavior of solutions to the corresponding limit problem given by and v x ∈ L ∞ loc ([0, ∞; L 2 ( )). (2.22) Proof As (2.21) consists of two actually decoupled scalar and linear drift-diffusion problems with Robin-type boundary conditions, global existence of a classical solution enjoying the claimed regularity properties is asserted by standard parabolic Schauder theory (Ladyzenskaja et al. 1968). To see its uniqueness within the indicated class, given two solutions (w, v) and w, v) both fulfilling the regularity requirements in (2.22), we write ϕ:=w − w and ψ:=v − v to see using (2.21) that and that ⎧ ⎨ (2.24) Therefore, by Young's inequality, for all t > 0 (2.25) and 1 2 Now for fixed t > 0, (2.22) asserts the existence of c 1 = c 1 (T ) > 0 and c 2 = c 2 (T ) > 0 such that v x (·, t) L 2 ( ) ≤ c 1 and w(·, t) L ∞ ( ) ≤ c 2 for all t ∈ (0, T ), whence in particular, using the Gagliardo-Nirenberg inequality and Young's inequality, we see that with some cx = c Thus, an integration using the initial conditions in (2.23) and (2.24) shows that ϕ 2 (·, t) + c 2 2 ψ 2 (·, t) = 0 for all t ∈ (0, T ) and therefore implies that ϕ ≡ ψ ≡ 0, for T > 0 was arbitrary.
Apart from that, thanks to the simple structure of the second equation in (2.21), one can readily achieve also some higher-order regularity information.
Lemma 2.7 If (1.8) holds, there exists C > 0 such that the solution of (2.21) satisfies Proof An application of standard parabolic Schauder theory ( Ladyzenskaja et al. (1968)) to the second subproblem contained in (2.21) provides θ 1 ∈ (0, 1) and (2.29) Therefore, when rewritten in the form w t = w x x + a(x, t)w x + b(x, t)w with a(x, t):= − v x and b(x, t):= − v x x enjoying suitable Hölder bounds according to (2.29), also the first equation in (2.21) becomes accessible to the same tool so as to allow for the existence of θ 2 ∈ (0, 1) and c 2 > 0 fulfilling Thanks to the latter, no nonzero boundary terms appear when we test (2.32) by − v x x and integrate by parts to see that for all t > 0.
(2.34) As (2.33) moreover enables us to invoke a Poincaré inequality to find c 1 > 0 fulfilling and that hence, upon integration, writing c 2 := v 2 x (·, 0), we have v 2 x (·, t) ≤ c 2 e −2c 1 t for all t > 0, (2.35) from which (2.30) immediately follows by definition of v and the fact that c 2 is finite due to (1.8). Secondly, a direct integration in (2.34) shows that where we once more use (2.32) to see that as a consequence of (2.35), In conclusion, (2.36) thus in particular entails that and thereby establishes (2.31).
Thanks to an evident mass conservation property of the second equation in (2.21), the latter in conjunction with one of the higher-order regularity features asserted by Lemma 2.7 already entails uniform stabilization of v toward one particular steady state: Lemma 2.9 Suppose that (1.8) holds. Then, there exists β > 0 such that the solution of (2.21) satisfies v(·, t) → β · exp in L ∞ ( ) as t → ∞. (2.37) Proof We let β:= v 0 e x dx and note that then β is positive by (1.8), and assuming (2.37) to be false, we could find (t k ) k∈N ⊂ (1, ∞) and c 1 > 0 such that t k → ∞ as k → ∞ and v(·, t k ) − β · exp L ∞ ( ) ≥ c 1 for all k ∈ N. (2.38) Now since (v(·, t)) t>1 is bounded in C 1 ( ) and hence relatively compact in C 0 ( ) by Lemma 2.7 and the Arzelà-Ascoli theorem, upon passing to a subsequence if necessary we may assume that we may rely on Lemma 2.8 to infer the existence of c 2 > 0 such that v ∞ (x) = c 2 for all x ∈ , (2.40) because (2.30) warrants that v kx → 0 in L 2 ( ) and k → ∞. By definition of v ∞ , in view of (2.39), this means that Since, on the other hand, a direct integration in (2.21) shows that this identifies c 2 according to c 2 = v 0 e x dx = β, so that (2.41) contradicts the hypothesis (2.38), thus altogether implying that actually (2.37) must have been valid.
Constituting the apparently most substantial part of this section, the following lemma turns the L 2 integrability property of v t contained in Lemma 2.8 into a first, though yet rather weak, information on stabilization in the first solution component.

Lemma 2.10
If (1.8) holds, then the solution of (2.21) has the property that x ∈ , t ≥ 0 and use the identity w x = e v ( w x + v x w), as thereby implied, to see that and that because of (2.21). Thanks to (2.43), an integration by parts in (2.44) shows that x for all t > 0.
(2.49) Therefore, namely, by utilizing Young's inequality, we can estimate (2.50) Since furthermore, again due to Young's inequality, and since another application of (2.49) shows that x for all t > 1, from (2.45) we all in all obtain that y(t):= e v(·,t) w 2 x (·, t) and h(t):= v 2 t (·, t), t ≥ 1, satisfy with c 7 := 1 c 6 e c 3 and c 8 := (c 2 +c 1 c 4 ) 2 c 6 e c 3 2 + 2. By straightforward integration, we thus infer that for all t > 2, because h(t) ≤ c 2 5 L for all t > 1 according to (2.48). Now since Lemma 2.8 ensures that this shows that which by definition of w yields (2.42) due to the fact that e v ≥ 1.
Again relying on Lemma 2.7, this can be seen to imply convergence to some equilibrium also in the first solution component:

Lemma 2.11
Assume (1.8). Then, there exists γ > 0 such that with β > 0 taken from Lemma 2.9, for the solution of (2.21), we have Proof Defining the number γ according to positive due to (1.8), we proceed as in Lemma 2.9 by firstly relying on Lemma 2.7 and the Arzelà-Ascoli theorem to see that (e −v(·,t) w(·, t)) t>1 is relatively compact in C 0 ( ), and by secondly identifying all corresponding ω-limits: Indeed, whenever (t k ) k∈N ⊂ (1, ∞) and ϕ ∈ C 0 ( ) are such that t k → ∞ and e −v(·,t k ) w(·, t k ) → ϕ in L ∞ ( ) as k → ∞, in view of Lemma 2.10 we have ϕ ≡ c 1 in with some c 1 ≥ 0, while Lemma 2.9 asserts that so that since independently we obtain from (2.21) that d dt w = 0 for all t > 0, we must have and hence c 1 = γ . As (t k ) k∈N was an arbitrary sequence having the indicated stabilization properties, by a standard reasoning in the style of that from Lemma 2.9, we readily infer (2.52).

Approaching (2.21) for Small Values of " in (1.7)
Our next goal consists in establishing a link between our solutions to (2.17) and those of (2.21) through an appropriate statement on convergence in the limit ε 0. The key ingredient toward this will be provided by the following outcome of a Moser-type iteration applied to the first equation in (2.17) on the basis of the taxic gradient estimate from (2.20).

Lemma 2.12
If (1.8) holds and T > 0, and if ε down (T ) > 0 is taken from Lemma 2.4, then one can find C(T ) > 0 such that the solution of (2.17) from Corollary 2.5 has the property that w ε (·, t) L ∞ ( ) ≤ C(T ) for all t ∈ (0, T ) and any ε ∈ (0, ε down (T )). (2.53) Proof For fixed T > 0 and nonnegative integers k, we let p k :=2 k and estimate by testing the first equation in (2.17) against w p k −1 ε for k ≥ 1 to see that due to Young's inequality, for all t > 0. Since herein εw ε ≤ 1 by (2.19), and since p k 2 ≤ p k − 1 ≤ p k for any such k, this shows that Now an application of Corollary 2.5 to q:=4 yields c 1 = c 1 (T ) > 0 such that v 4 εx ≤ c 1 for all t ∈ (0, T ) and any ε ∈ (0, ε down (T )), whence on the right-hand side of (2.54), we can use the Cauchy-Schwarz inequality to estimate Here, by means of the Gagliardo-Nirenberg inequality, we can find c 3 = c 3 (T ) > 0 such that due to our definition of (M j ) j≥0 and Young's inequality, and therefore entails that and ε ∈ (0, ε down (T )), thus implying (2.53) also in this case.
Indeed, thanks to a known result from scalar parabolic theory, the latter, once more in conjunction with (2.20), ensures bounds for both solution components in (2.17) in appropriate Hölder spaces.
In quite a straightforward manner, the compactness features gathered above, along with the uniqueness statement contained in Lemma 2.6, enable us to take ε 0 in (2.17) with the desired result.

Lemma 2.15
Suppose that (1.8) is valid, and for T > 0 let ε down (T ) > 0 be as given by Lemma 2.4. Then the solutions of (2.17) gained in Corollary 2.5 have the properties that where (w, v) denotes the unique classical solution of (2.21) satisfying (2.22).
Proof In view of the uniqueness statement in Lemma 2.6, according to a standard argument, it is sufficient to make sure that with ε down (T ) > 0 taken from Lemma 2.4, each sequence (ε j ) j∈N ⊂ (0, ε down (T )) such that ε j 0 as j → ∞ contains a subsequence (ε j k ) k∈N such that (2.62)-(2.64) hold as ε = ε j k 0, with some classical solution (w, v) of (2.21) fulfilling (2.22). To verify this, given any such (ε j ) j∈N , we may rely on the bounds provided by Lemmas 2.13, 2.14 and Corollary 2.5 to see upon a straightforward extraction procedure involving the Arzelà-Ascoli theorem that in fact it is possible to find a subsequence (ε j k ) k∈N and functions w and v for which (2.62)-(2.64) hold as ε = ε j k 0. Since from (2.62) and (2.63), it is an evident consequence of the inequality κ > 1 that the validity of (2.21) for these functions results in taking ε = ε j k 0 in each of the summands in (2.17) in the classical pointwise sense, we furthermore conclude that (w, v) indeed forms a classical solution of (2.21), having the additional feature (2.22) due to (2.64), whereby the proof becomes complete.

Proof of Theorem 1.1
We are now prepared to derive our main result on preferred downstream migration in the presence of suitably small population densities by simply combining Lemma 2.15 with Lemma 2.11. proof of Theorem 1.1 Given δ > 0, by means of Lemma 2.11, we first choose T (δ) > 0 large enough such that the solution of (2.21) satisfies for all x ∈ and t ≥ T (δ). (2.65) We then fix any T > T (δ) and take ε down (T ) > 0 as accordingly provided by Lemma 2.4, and noting that then classical solvability of (1.7) in × (0, T ) within the desired class is asserted by Lemma 2.4, we may rely on Lemma 2.15 to find some ε 0 (δ, T ) ∈ (0, ε down (T )) such that for the correspondingly obtained classical solution, we have for all t ∈ [0, T ] whenever ε ∈ (0, ε 0 (δ, T )).
(2.66) and since another application of the maximum principle warrants uniqueness of classical solutions from X × X to these actually decoupled scalar transport-diffusion problems, this already implies the convergence property (z ε , v ε ) → (z, v) along the entire net (0, ε up ) ε 0, with (z, v) solving (3.23). The claim of the lemma is an evident by-product thereof.
As a consequence, we infer that the behavior of solutions to (1.11) is essentially determined by that of solutions to (3.7) in the sense specified in our main statement on dominance of upstream migration: proof of Theorem 1.2 Proceeding as in the proof of Theorem 1.1, for fixed δ > 0, we first invoke Lemmas 3.3 and 3.2 to pick T (δ) > 0 such that for the solution of the limit problem (3.7), abbreviating β:= κ 2 4 + π 2 L 2 and taking α > 0 from (3.11), we have which in particular entails that for all x ∈ and t ≥ T (δ). (3.24) Then, letting T > T (δ) be given, we can thereafter rely on Lemma 3.6 in choosing ε 0 (δ, T ) ∈ (0, ε up ) small enough such that e βt z ε (x, t) − e βt z(x, t) ≤ δ 2 for all x ∈ and t ∈ [0, T ], (3.25) which when combined with (3.24) yields (1.13) with γ :=α.

Numerical Investigations
We would like to carry out some numerical investigations that can illustrate the essence of Theorems 1.1 and 1.2. For that purpose we replace the model (1.5) by a discrete approximation. We refer to "Appendix A" for details. For the parameters in (1.6), we set κ = 2 = λ. Then, we get The function g(u) is introduced to rewrite (1.5) in a form convenient for the discretization presented in "Appendix A."

Case 1: Sparsely Distributed Cell Population
As initial data w 0 (x) and v 0 (x) to be used in (1.7), the following functions are specified  figure) are close to each other but rather far from the steady state w. Bottom: T = 0.5. The solution w of the limit problem is close to the steady state w; however, there is a certain discrepancy between w ε and the steady state w, which is natural to be attributed to the relatively large ε > 0 (4.1) We recall that the analytical steady-state solutions v(x) and w(x) of the limit problem (2.21) are found to be, see (2.37) and (2.52), (4.2) We consider a numerical grid with N = 100 grid cells on the spatial domain [0, 1]. We are interested in an illustration of the result of Theorem 1.1. We set ε = 0.5. We see that for this choice of ε > 0, the discrepancy between w ε and the steady state w seen in Fig. 4 has essentially been eliminated We compute numerical approximations of the solution (u ε , v ε ) of the initial boundary value problem (1.7) at a specified time T . We visualize w ε (x, T ) = ε −1 u ε (x, T ) and v ε (x, T ), which naturally allow us to compare with numerical approximations to w(x, T ) and v(x, T ), i.e., the solution of the limit problem (2.21), as well as the analytically computed steady states w(x) and v(x). The results are shown in Fig. 4. Focusing on w ε , the main observation is that for this choice of ε, there is a visible difference between w ε and the steady state w that remains. According to Theorem 1.1, this difference is controlled by ε. To test this, we reduce it by setting ε = 0.1. The results are shown in Fig. 5 and clearly illustrate the essence of Theorem 1.1 as expressed by the precise estimate in (1.10): For a sufficiently small ε > 0, we can ensure that the solution w ε = u ε /ε is as close to the steady state w as we want, subject to the condition that t is sufficiently large. . We see that the discrepancy between e βt z ε and the analytical steady state αe − κ 2 x sin π x is controlled by ε

Case 2: Densely Distributed Cell Population
We are interested in an illustration of the result of Theorem 1.2. For that purpose, we compute numerical approximations of the solution (u ε , v ε ) of the initial boundary value problem (1.11) at a specified time T = 0.1, but now with initial data u ε (x, t = 0) = 1 − εw 0 . First, we set ε = 0.5, i.e., our initial data are close to 1 apart from a local region in the center of the domain. We compute the corresponding solution u ε by using the numerical scheme for (1.5) specified in "Appendix A" and then compute the corresponding z ε using the relation u ε = 1 − εz ε . This allows us to compare e βt z ε (x, T ) and the steady-state expression αe − κ 2 x sin π x where α is given by (3.11), β = κ 2 4 + π 2 , and we use that L = 1. We want to show that we can control this discrepancy by the choice of ε, consistent with estimates (3.24) and (3.25). In Fig. 6, results are shown for ε = 0.5 (left) and ε = 0.1 (right). The numerical simulations indicate in accordance with Theorem 1.2 that by choosing ε sufficiently small and time T sufficiently large, we ensure that the solution u ε is as close as we want to an analytical steady-state solution that is skewed in the upstream direction.

Appendix: A Discrete Scheme for the Model (1.5)
We may write model (1.5) in the following form: (4.3) We introduce the notation g(u) = (4.7) where D + a j = a j+1 −a j x and h n j+1/2 := u n j ψ(u n j+1 ) when ([v − g(u)] n j+1 − [v − g(u)] n j ) ≥ 0 u n j+1 ψ(u n j ) when ([v − g(u)] n j+1 − [v − g(u)] n j ) < 0.  A similar discrete scheme can be defined for the limit problem (2.21) which allows us to compute discrete approximations {w n j } and {v n j }.