Traveling waves in a coarse‐grained model of volume‐filling cell invasion: Simulations and comparisons

Many reaction–diffusion models produce traveling wave solutions that can be interpreted as waves of invasion in biological scenarios such as wound healing or tumor growth. These partial differential equation models have since been adapted to describe the interactions between cells and extracellular matrix (ECM), using a variety of different underlying assumptions. In this work, we derive a system of reaction–diffusion equations, with cross‐species density‐dependent diffusion, by coarse‐graining an agent‐based, volume‐filling model of cell invasion into ECM. We study the resulting traveling wave solutions both numerically and analytically across various parameter regimes. Subsequently, we perform a systematic comparison between the behaviors observed in this model and those predicted by simpler models in the literature that do not take into account volume‐filling effects in the same way. Our study justifies the use of some of these simpler, more analytically tractable models in reproducing the qualitative properties of the solutions in some parameter regimes, but it also reveals some interesting properties arising from the introduction of cell and ECM volume‐filling effects, where standard model simplifications might not be appropriate.


INTRODUCTION
Cell invasion is central to a variety of biological phenomena, playing a key role in morphogenesis, tumor growth, and tissue engineering.Many different mathematical approaches have been used to model cell invasion processes, from agent-based models describing the processes underlying invasion at a single-cell level, to partial differential equation (PDE) models that provide a cell population-level description of invasion in terms of cell density dynamics. 1While some of these PDE models have been formulated through adaptations of classical models for invasion processes in other biological contexts, many are derived by coarse-graining a cell-level description to produce PDEs that offer the corresponding population-level description.However, there remain a number of unanswered questions regarding the specific choice of model for a given application, including how varying assumptions in the agent-based model give rise to different PDE models and to what extent these differences impact the cell population-level description. 2n this work, the impact of various modeling assumptions at the single-cell level is compared by investigating the qualitative and quantitative properties of the solutions of the resulting PDE models.In particular, since in many real-life instances of cell invasion, the cells have to invade through extracellular matrix (ECM) [3][4][5] -that is, the network of proteins and other molecules that impact collective cell invasion by reducing space available for cells to migrate into-the focus is on models of cell invasion into ECM.
The classic example of a model describing invasion of a single population is the Fisher-Kolmogorov-Pietrovskii-Piskunov (FKPP) model, which was first proposed in the context of the spread of an advantageous gene. 6,7This model has seen a broad spectrum of applications in the natural sciences: most notably in cell biology 8,9 and ecology, 10,11 where traveling wave solutions are representative of invasion phenomena. 12,13A model of cell invasion through ECM is presented in El-Hachem et al. 14 This consists of a system of two coupled PDEs, with a nonlinear cross-species density-dependent diffusion term and logistic growth, whereby proliferation of the cell population is limited by the presence of cells and ECM.A similar model is considered in Colson et al., 15 where proliferation depends only on the presence of cells.An obvious question to ask is how the predictions of such models may be affected by a consistent description of the role of volume-filling effects (i.e., cells and ECM take up some given volume, preventing cell invasion) across both proliferative and diffusive mechanisms of cell dynamics.In particular, the impact of crowding on cell motility is generally modeled at a population-level by the introduction of a density-dependent diffusion term, such as in Refs.[16-19].However, these models provide a phenomenological description of the impact of crowding, rather than considering how interactions at the individual cell-level directly impact motility at the population-level.
This study aims to extend and apply the work in Simpson et al. 20,21 to develop an agent-based model for cell invasion into ECM, taking into account volume-filling effects, where both cell motility and proliferation are impacted by the presence of other cells or ECM components.We make the simplifying assumption that space is the only factor limiting cell invasion, whereas other models in the literature 22,23 assume various other factors, such as nutrient-limited growth. 24By coarsegraining this model, a limiting PDE description is formally derived and explored both analytically and numerically, making it possible to carry out a systematic comparison between the populationlevel behaviors observed in this model and those predicted by simpler models.In particular, we compare how the population-level behaviors predicted by this model relate to existing models built on different constitutive assumptions, such as the FKPP model, or the simpler models presented in Refs.[14, 15].Each of these simpler models can be recovered from the model presented in this work by neglecting specific terms, such as those capturing volume-filling effects.

MATHEMATICAL MODEL AND PRELIMINARY RESULTS
We begin by developing a simple one-dimensional, on-lattice, agent-based model of cell invasion into ECM that incorporates both cell motility and proliferation, and degradation of ECM, in the presence of volume-filling effects.We then coarse-grain this model to formally derive a corresponding PDE model that comprises a system of coupled PDEs for the densities of cells and ECM. 20,25

Agent-based model
In the simplified setting of this model, cells are represented as discrete agents that can proliferate and move on a one-dimensional uniform lattice, which constitutes the spatial domain and can also degrade the surrounding ECM, which is regarded as being composed of discrete constitutive elements.The novel aspect of this model is the introduction of volume-filling effects, similar to the model described in Morris et al., 26 which uses the methods described in Painter and Hillen, 27 but extended to multiple populations. 21et the number of cells and ECM elements on lattice site  = 1, 2, … of width Δ at time t ∈ ℝ + of realization  = 1, 2, … ,  of the model be denoted, respectively, by    ( t) and    ( t).We assume that ECM elements have constant density and are chosen to occupy a volume equal to that of a cell, such that at most  cells or ECM elements can occupy each lattice site.
The dynamics of the cells are governed by two mechanisms: proliferation, in which a cell places a daughter cell into the same lattice site it occupies; and motility, whereby cells can move to one of their two adjacent lattice sites.Moreover, ECM elements can be degraded by cells in the same lattice site as them.To incorporate volume-filling effects into the model, we prescribe that each lattice site has a maximum occupancy level , 28 and we assume that (A1) if a cell attempts a move to a neighboring lattice site, then the probability that the move is successful decreases linearly with the occupancy level of the target site, such that the probability of a successful move to a target site with occupancy level  is zero; (A2) if a cell attempts to proliferate, then the probability of success decreases linearly with the occupancy level of the site where the cell is located, such that the probability of a successful proliferation event in a site with occupancy level  is zero.

Probability of cell movement
A cell attempts a movement in a time step  with probability  m ∈ [0, 1], and the attempted movement from lattice site  to either of the neighboring lattice sites  ± 1 occurs with equal probability 1∕2.Using Assumption 2.1, we can define the probability of movement to the left,  m  − ( t), or right,  m  + ( t), during the time interval [ t, t + ) of realization , as (

Probability of cell proliferation
Note that, the initial distributions of cells and ECM elements must be such that at most  cells or ECM elements can occupy each lattice site to ensure the probabilities  m  ± ( t),  p   ( t) ≥ 0 are well-defined.Under the assumption that the initial distributions of cells and ECM elements satisfy 0 ≤    (0) +    (0) ≤  for all  = 1, 2, … ,  and  = 1, 2, …, the definitions for the probabilities of cell movement and proliferation given by Equations ( 1) and (2) ensure that 0 ≤    ( t) +    ( t) ≤  for all  = 1, 2, … ,  and  = 1, 2, … for any t ∈ ℝ + . (3)

Probability of ECM degradation
During the time interval [ t, t + ) of realization , an element of ECM in lattice site  is degraded by a cell on the same lattice site with probability  d ∈ [0, 1], such that the degradation per unit element of ECM,  d   ( t), is

Corresponding coarse-grained model
In order to derive a coarse-grained description of the agent-based model, we introduce the average occupancy of lattice site  at time t by cells and ECM elements over  realizations of the model, denoted, respectively, by

Coarse-grained model of cell dynamics
We proceed to derive a coarse-grained model by considering how the average occupancy in lattice site  changes during the time interval [ t, t + ): Note that, in writing down Equation (4) we have used probabilistic approximations of the meanfield type which are frequently used for the coarse-graining of agent-based models and involve assuming independence of lattice sites (see, e.g., Penington et al. 29 ).Rearranging Equation ( 4) and dividing both sides by  yields We now divide both sides of Equation ( 5) by length scale Δ, perform a Taylor expansion, and take limits as Δ,  → 0 to obtain a description of the cell density dynamics in terms of the variables ũ( x, t) and m( x, t), that are the continuum counterparts of ⟨  ( t)⟩∕Δ and ⟨  ( t)⟩∕(Δ) that represent, respectively, the number density of cells and the density of ECM at position x ∈ ℝ and time t ∈ (0, ∞).The factor  represents the number of cells equivalent to a unit mass of ECM and is introduced as a conversion factor between the density of ECM, as defined by mass of ECM per unit volume, and the number density of ECM elements, given by  m( x, t).Under the assumptions we obtain the following PDE for the cell density ũ( x, t): where x ∈ ℝ and t ∈ (0, ∞).Note that the first term on the right-hand side of Equation (7)  describes the movement of cells down gradients in cell density, with movement prevented by the presence of surrounding cells and ECM, as expected by the introduction of volume-filling effects.
The second term models the motion of the cells down the "total density gradient" of cells and ECM, ũ + μ m.The third term captures cell proliferation, which is also impacted by volume-filling effects.From Equation (7), it is clear that the parameter D ≥ 0, which is defined via Equation ( 6), can be regarded as the diffusion coefficient of the cells in the absence of ECM, while the parameters r ≥ 0 and K > 0, which are also defined via Equation (6), are the intrinsic growth rate of the cell population, and the density corresponding to the maximum occupancy level (i.e., the carrying capacity), respectively.

Coarse-grained model of ECM dynamics
Probabilistic approximations similar to those underlying Equation (4) give the following conservation equation for the evolution of ECM elements in lattice site  during the time interval [ t, t + ): Rearranging Equation ( 8), dividing by Δ and  and taking limits as Δ,  → 0, under the assumption we formally obtain the following differential equation for ECM density m( x, t): where x ∈ ℝ and t ∈ (0, ∞).Here, the parameter λ ≥ 0 defined via Equation ( 9) is the per cell degradation rate of ECM.We observe that when there is no ECM degradation (i.e., if λ = 0) and ECM is uniformly distributed at t = 0 (i.e. if m( x, 0) ≡ m0 where m0 ∈ ℝ + with 0 ≤ m0 ≤ K), the mathematical model defined via Equations ( 7) and (10) simplifies to the following FKPP model of cell dynamics 6 : where

Nondimensional coarse-grained model
The mathematical model defined via Equations ( 7) and ( 10) can be nondimensionalized by the introduction of the following nondimensional variables and written as where  ∈ ℝ and  ∈ (0, ∞).Here, the only remaining parameter is  = λ K∕ r ≥ 0, which is interpreted as the rescaled ECM degradation rate.We complement the model defined via Equations ( 12) and ( 13) with no flux boundary conditions for Equation ( 12): and , ∕ → 0 as  → ∞.We also have the following initial conditions: We note that by assuming at the single-cell level that both the presence of cells and ECM elements impair the movement and proliferation of the cells, the resulting population-level description for cell density evolution in Equation ( 12) exhibits a number of differences to similar models without volume-filling effects.For example, the model studied by El Hachem et al. 14 does not consider volume-filling of cells to impair cell movement, and therefore contains one less flux term, namely that accounting for movement of cells down the "total density gradient."This model can be recovered from Equation ( 12) by employing different underlying assumptions such that the probability of movement depends on the average available space (where space is only filled by ECM) between the target lattice site and the lattice site the cell occupies at time t.

Numerical exploration of possible traveling wave solutions
We are interested in the possible constant profile, constant speed traveling wave solutions displayed by the model defined via Equations ( 12) and ( 13).As such, we first explore the range of possible behaviors numerically.We report on the results of numerical simulations carried out for the model posed on the spatial domain (0, ), with  > 0 sufficiently large so that the no flux boundary condition ( 14) at  =  does not interact with the traveling wave.The simulations were subject to the following initial conditions: F I G U R E 1 Numerical solutions of Equations ( 12) and ( 13) subject to the initial conditions ( 16) and ( 17

𝑚(𝑥
where 0 <  ≪  represents the width of the initially invaded region at  = 0 and  0 ∈ [0, 1) corresponds to the uninvaded density of ECM ahead of the cells.We note that, by design, the model ( 12) and ( 13) does not permit traveling waves when there are initial conditions with compactly supported cell density and  0 = 1.This is because cells require space ahead of the wave in order to invade; in any regions initially devoid of cells, the ECM cannot be degraded to allow cells to invade.As such, we proceed by considering  0 ∈ [0, 1).Further specifics of the parameter values and the numerical methods used in this paper can be found in Appendix B.
As shown in Figure 1, when  0 ∈ [0, 1) the solutions to Equations ( 12) and ( 13) subject to the initial conditions ( 16) and ( 17) converge to traveling waves whereby the cell density, , decreases monotonically from one to zero and the ECM density, , increases monotonically from zero to  0 .The numerical results in Figure 1 also indicate that the speed of the traveling waves changes as the values of the parameters  and  0 are changed.This is illustrated in more detail in Figure 2 that also shows that (in agreement with the analytical results presented in Section 3) when  0 ∈ (0, 1): if  → 0 + then the speed of the traveling waves converges to  = 2 (1 −  0 ); whereas if  → ∞ then the speed of the traveling waves converges to  = 2.
We also note that when  0 = 0, the solutions to Equation ( 13) subject to the initial condition ( 17) are such that (, ) ≡ 0 for all  ≥ 0 and thus the model simplifies to the FKPP model (11) with Consistent with this, numerical simulations indicate that when  0 = 0, the cell density  converges to a traveling wave that decreases monotonically from one to zero (results not shown), and travels with speed  = 2 (i.e., the minimal speed of traveling wave solutions to the FKPP model ( 18)) (see Figure 2).

F I G U R E 2
The relationship between the numerically estimated speed (solid lines) of traveling wave solutions of Equations ( 12) and ( 13) subject to the initial conditions ( 16) and ( 17).The dashed lines in the plot on the left highlight the value of 2(1 −  0 ).The numerically estimated traveling wave speed is obtained by tracing the point () such that ((), ) = 0.1.Further specifics of the parameter values and the numerical methods used can be found in Appendix B.
The numerical results summarized in Figure 2 for  0 ∈ [0, 1) show similar behaviors to that in El-Hachem et al., 14 where no volume-filling effects of cells prevent cell movement, while a marked difference is observed for the case  0 = 1, as discussed in Appendix A.
The steady states of the system ( 29) and ( 31) with boundary conditions ( 32) and ( 33) are given by  1 = (1, 0, 0) and  2 = (0, 0,  0 ).Traveling wave analysis based on standard linear stability techniques (i.e., standard traveling wave analysis) 30 seeks trajectories in the phase space that connect  1 at  = −∞ to  2 at  = ∞. 31,32The eigenvalues of the linearized system at (, , ) = (1, 0, 0) are which implies that (1,0,0) is a three-dimensional, hyperbolic, unstable saddle point, 33 which has eigenvectors given by F I G U R E 3 Phase plane plot of the system of ordinary differential equations ( 29) and ( 30), for different traveling wave speeds, , demonstrating the change from a stable spiral to a stable node as the traveling wave speed exceeds  min .The corresponding unstable eigenvector given by Equation ( 36) and stable eigenvector given by Equation ( 39) are overlaid in the lower plots. 30Further specifics of the parameter values and the numerical methods used can be found in Appendix B.
, = The eigenvalues of the linearized system at (, , ) = (0, 0,  0 ) are with corresponding eigenvectors which implies that (0, 0,  0 ) is a stable, nonhyperbolic fixed point 34 (see Appendix D for the full derivation).In all cases, we use index 2 to refer to the positive of the two choices and 3 for the negative.When  2 − 4(1 −  0 ) 2 ≤ 0, the steady-state (0, 0,  0 ) is a stable spiral as the eigenvalues have nonzero imaginary parts; however, when  2 − 4(1 −  0 ) 2 ≥ 0, the steady state is a stable node.In the case that the state (0, 0,  0 ) is a stable spiral,  oscillates around this point on its approach and can therefore take negative values (see Figure 3).However, when (0, 0,  0 ) is a stable node, there must exist a trajectory from (1,0,0) to (0, 0,  0 ) contained entirely in the region of phase space defined by  ≥ 0,  ≤ 0, and  ≥ 0, which ensures nonnegativity of  and , as required to be biologically consistent.This demonstrates the existence of a minimum wave speed,  min = 2(1 −  0 ), such that the dependent variables,  and , remain nonnegative for all time.It is important to note that  min is a lower bound on the traveling wave speed, which is only actually attained for this system when the rescaled ECM degradation rate is sufficiently small, that is,  → 0 + (see Section 3.1).This is clearly shown in Figure 2, which also demonstrates that decreasing  0 results in an increase in the traveling wave speed.
Traveling wave analysis has also been performed on a PDE model for melanoma invasion into the skin, 35 where volume-filling effects of cells are not considered to impact cell movement, 14 as described earlier.Since traveling wave analysis is always performed on the linearized system, it follows that the additional term describing cell movement prevented by the presence of other cells is lost from Equation (30) during linearization and the minimum traveling wave speed is the same as that derived in El-Hachem et al. 14 :  min = 2(1 −  0 ).Another minimal model for tumor growth was proposed in Colson et al., 15 where the volume-filling effects of cells were not accounted for in describing cell movement or cell proliferation.Both of these models have the same equation for ECM density as Equation ( 13), and the models in Refs.[14, 15] have the same flux terms in the equation for cell density evolution, but the model in Colson et al. 15 has one less reaction term since proliferation is unimpeded by the local ECM density.As a result of the fact that all volume-filling effects are encoded in nonlinear terms, changes to the flux terms alone (within this suite of models) have no effect on the predicted minimum traveling wave speed, as they are all identical after linearization.However, alterations to the net proliferation terms do significantly impact the minimum traveling wave speeds predicted by standard traveling wave analysis.Further information regarding these models and their differences can be found in Appendix C.
As previously described, we are particularly interested in investigating the dependence of traveling wave solutions on the parameters , the rescaled ECM degradation rate, and  0 , the density of ECM far ahead of the wave.Having now determined that the minimum traveling wave speed decreases linearly as  0 increases, we now aim to explore the relationship between the numerically estimated traveling wave speed and .
Since the traveling wave speed depends on , standard perturbation techniques are difficult to apply to the traveling wave equations (19) and (20).As a result, we examine Figure 2 for clues as to how to proceed.We immediately see that for sufficiently small  it appears that the numerically estimated traveling wave speed is independent of  and matches the speed predicted by standard traveling wave analysis.It can also be seen from the contour plot in Figure 2 that for large values of , the speed converges for all values of  0 ∈ [0, 1).As such, we now investigate the asymptotic limits corresponding to slow and fast rescaled ECM degradation rates,  → 0 + and  → ∞, respectively.

Formal asymptotic analysis for 𝝀 → 𝟎 +
Using Equation (28), it is clear that for  ∈ (, ∞) (see Figure 4 or Figures E1 and E2 for the traveling wave profiles).
In the asymptotic regime  → 0 + , substituting Equation (40) into Equation ( 24) and using the fact that, since 0 ≤ () < 1 for  ∈ (, ∞) and ()∕ ≈ − () for  ∈ (, ∞) (cf. the ansatz F I G U R E 4 Numerical solutions of Equations ( 12) and ( 13) subject to the initial conditions ( 16) and ( 17 given by Equation ( 27)), the following asymptotic relation holds: for  ∈ (, ∞), we find for  ∈ (, ∞).Equation ( 42) is equivalent to the FKPP model (11) in traveling wave coordinates with D = r = K = 1 −  0 , with ĉmin = 2(1 −  0 ), as predicted earlier.An excellent match around the leading edge of the traveling wavefront between the traveling wave solution to the FKPP model (11) (Equation ( 43) in traveling wave coordinates) and Equations ( 12) and ( 13) for low values of the rescaled ECM degradation rate can be seen in the plot on the left in Figure 6 (see also Figure 4 or Figures E1 and E2).We now consider the region  ∈ (−∞, ) by rescaling Equations ( 19) and (20) using the new variable  =  for  ∈ (−∞, ].The system of Equations ( 19) and (20) becomes F I G U R E 5 Numerical solutions of Equations ( 12) and ( 13) subject to the initial conditions ( 16) and ( 17), for  0 = 0.2 on the left and  0 = 0.8 on the right with rescaled ECM degradation rate  = 10 −3 translated into the traveling wave coordinate, .Solid lines represent the cell and ECM densities from numerical simulations in purple and orange, respectively.The FKPP solution (42) in traveling wave coordinates is plotted as a dotted blue line.The solution () =  0 is plotted in dotted black, and the analytical solutions given by Equations ( 48) and ( 49) are plotted in dashed blue and black lines, respectively.Further specifics of the parameter values and the numerical methods used can be found in Appendix B.

F I G U R E 6
Left: plot of the cell density, , obtained through numerical simulations of Equations ( 12) and ( 13) subject to the initial conditions ( 16) and ( 17) (solid lines) for small values of , and numerical simulations of the FKPP model (11) with rescaled coefficients D = r = K = 1 −  0 (dashed black line) with  = 100 and  0 = 0.6.Right: plot of the cell density, , obtained through numerical simulations of Equations ( 12) and ( 13) subject to the initial conditions ( 16) and ( 17) (solid lines) for large values of , and numerical simulations of the FKPP model (18) (dashed black line) in the plot on the right for  = 50 and  0 = 0.4.Qualitatively, the same behavior is observed for all  0 ∈ [0, 1).Further specifics of the parameter values and the numerical methods used for the simulations can be found in Appendix B.
In the traveling wave coordinate, , the solution for the wave profile of the ECM given by Equation ( 46) is for  ∈ (−∞, ], and () =  0 for  ∈ (, ∞), as given by Equation ( 40).An excellent agreement between these analytical solutions and the numerical results can be observed in Figure 5.
Similar models, such as those described at the end of Section 3 that do not have volume-filling effects taken into account, demonstrate qualitatively similar behavior.In all of these models, at very low rescaled ECM degradation rates we observe convergence of the solutions to those of the FKPP model with rescaled parameters.For models with the same cell proliferation term as in Equation ( 12), the rescaled parameters are the same and the convergence has qualitatively similar behavior, as displayed in the plot on the left in Figure 6.As a result, in the limit of very small rescaled ECM degradation rates,  → 0 + , the model ( 12) and ( 13) can be simplified to that presented in El-Hachem et al., 14 which neglects the volume-filling effects of cells upon cell movement.This model can, in turn, be well approximated by the FKPP model (11) with rescaled parameters D = r = K = 1 −  0 .This result is consistent with predictions from standard traveling wave analysis.However, for the model presented in Colson et al., 15 the parameters of the rescaled FKPP model to which the model converges are, instead, D = 1 −  0 and r = K = 1, that entails a higher cell carrying capacity density since proliferation is not impacted by the surrounding ECM.(See Appendix C for a more detailed comparison.)As such, the model ( 12) and ( 13) is poorly approximated using models, such as that in Colson et al., 15 with different underlying assumptions for cell proliferation.These differences highlight the importance of fully laying out all of the model assumptions at the single-cell level before deriving the PDE model, so that the populationlevel model fully captures behaviors associated with the underlying cell-level assumptions, in all parameter regimes.

Formal asymptotic analysis for 𝝀 → ∞
In the case of very large rates of ECM degradation, by considering the semiexplicit solution for  in terms of  given by Equation ( 28), we see that for  ∈ (, ∞) (see Figure 7 or Figure E3 for the traveling wave profiles).In the asymptotic regime  → ∞, substituting Equation (50) into Equation ( 24) and using the fact that, since 0 ≤ () < 1 for  ∈ (, ∞) and ()∕ ≈ − () for  ∈ (, ∞) (cf. the ansatz given by Equation ( 27)), the Numerical solutions of Equations ( 12) and ( 13) subject to the initial conditions ( 16) and ( 17 following asymptotic relation holds: for  ∈ (, ∞), we find for  ∈ (, ∞).Hence, when  → ∞ we expect () at the leading edge of the traveling front to behave, to a first approximation, as the solution to the FKPP equation ( 18) in traveling wave coordinates subject to the boundary condition (22), for which  min = 2.This result can also be observed numerically in the plot on the right of Figure 6.The same behavior is observed in similar models without volume-filling effects, 14,15 demonstrating that the model ( 12) and ( 13) can be approximated, to an extent, with any of these simpler models in the parameter regime  → ∞, as growth and diffusion are unrestricted by the ECM within a neighborhood of the traveling wavefront.

DISCUSSION AND CONCLUSIONS
In this paper, a model for cell invasion into the surrounding ECM has been studied by considering primarily its traveling wave solutions.In this model, derived from the first principles from an agent-based model describing cell-level behaviors, cells evolve under the action of diffusion and proliferation, that is coupled to degradation of the surrounding ECM.As a result of volume-fillling effects, cells require space ahead of the wavefront in order to invade the domain.Numerical solutions of the PDE model ( 12) and ( 13) demonstrate a complex relationship between the traveling wave speed, , the density of ECM far ahead of the wave of cells,  0 , and the rescaled ECM degradation rate, .Partial relationships between these parameters in asymptotic regimes of interest have been established, including that  → 2(1 −  0 ) as  → 0 + and that TA B L E 1 Description of the volume-filling effects of cells and ECM considered by the models compared in this study.Equations ( 12) and ( 13)
→ 2 − as  → ∞.A good agreement with the FKPP model (11) has been demonstrated in the case where  → ∞, and we showed that the impacts of introducing volume-filling effects of cells to reduce cell movement (in comparison to the model in El-Hachem et al. 14 ) are minimal.As such, the FKPP model ( 11) provides a suitable model simplification to reproduce the qualitative behaviors of the fully dimensional system in the case of a large ECM degradation rate, λ, compared to the proliferation rate, r.Since  = λ K∕ r, the results equivalently suggest that as K → ∞, the system can be well modeled by the FKPP model (11).This describes a model where volumefilling effects are negligible, and thus the speed of the invasion front is given by  min = 2.For  → 0 + , which is representative of very large proliferation rates compared to the rescaled ECM degradation rates, or extremely small carrying capacities, the system can be studied by considering the simplification to a rescaled FKPP model (42).In this case, traveling waves are observed for  ∈ [0, 1), but the speed of the invasion front is now given by  min = 2(1 −  0 ).Converting back to dimensional variables, as with the FKPP model (11), the analytically predicted traveling wave speed increases with the cell proliferation rate, but with a more complicated relationship for the regions of parameter space corresponding to where the relationship between the traveling wave speed and rescaled ECM degradation rate is not yet well established.It is likely this complicated relationship indicates that the system exhibits changes between pulled, pushed, and semipushed waves due to the nonlinear cross-species dynamics that vary in strength for different parameter values. 36This could be investigated further by examining the ratio between the traveling wave speeds for different parameter values.
It is also clear that qualitatively similar results are observed between this new model with volume-filling, and previously studied models outside this framework, as described by Table 1, in all cases where  0 ∈ [0, 1).Therefore, it could be said that the model originally proposed in Browning et al. 35 provides a good model simplification for any case where  0 ∈ [0, 1).In the case where  0 = 1, the region that is initially uninvaded by cells is full of ECM, such that proliferation and movement of cells into this region is entirely prevented.This result provides the starkest difference between the model studied in this paper and those previously studied elsewhere. 14,15It is observed that in the case of compactly supported initial cell density, cell invasion cannot occur into the region where (, 0) = 1, (, 0) = 0, and thus pinning occurs and traveling waves cannot form. 37It is biologically reasonable to assume that an invading cell population might have zero density far ahead of the invading front.However it is important to note that the model considered here is a very simplistic model for cell invasion into ECM, and if further biological complications, such as the secretion of matrix metalloproteinases (MMPs) by cells to degrade and remodel ECM, were introduced then these phenomenological results would no longer be observed. 38This is because we could reasonably assume MMPs could still diffuse into regions occupied entirely by ECM, and then degrade it.
The overall conclusion of our study is that there exist simpler models for cell invasion into ECM such as Refs.[6, 14, 15] that are defined by similar guiding principles and can be used to reproduce the qualitative behaviors of the traveling waves observed in the model presented in this work.Analysis of these systems confirms that the qualitative model predictions are conserved, and therefore the simpler models can be used in future studies to reduce computational complexity and make the resulting PDE model more analytically tractable.The disadvantage of this conclusion, however, is that in order to use these models to infer parameters from data, extra steps would be required to validate whether the correct model has been selected.0][41] Our results reveal that the reaction term significantly impacts the traveling wave speed for small and intermediate values of  and thus it could be used to inform model development, by defining the reaction term by considering whether space or nutrients are the limiting factor for cell invasion into ECM; and model selection, by comparing the expected wave speeds to the data.
There are a variety of possible extensions to the work presented in this paper.The underlying on-lattice agent-based model of cell movement involves a number of simplifying assumptions, such as that cells can only degrade ECM agents in the same lattice site.By varying these assumptions, there would be the possibility to expand the biological applicability of the study to determine under which regimes the resulting models can also be approximated by simpler seminal models of cell invasion.Different proliferation terms, as well as terms to account for ECM evolution in more detail could be included, such as ECM remodeling by cells, or elastic deformation. 42Beyond this, another clear extension of this work would be to introduce further spatial dimensions, or different geometries, that are particularly interesting for studying cancer cell invasion, and to investigate the stability of the traveling wave solutions for the different possible models.For the case  → 0 + , there is an opportunity to apply boundary layer theory and asymptotic analysis to arrive at an expression for the full traveling wave profile at long times.It would also be of particular interest to arrive at some functional form for the traveling wave speed, (,  0 ), for all possible parameter values and to define the critical value of   , depending on  0 (see Figure 2), whereby for  <   the minimum traveling wave speed observed numerically matches that predicted by standard traveling wave analysis  =  min .The critical value,   , might be found by establishing the basins of attraction for each steady-state and seeking parameter regimes where the dynamics follow different paths.If possible, this knowledge could then further aid an investigation using perturbation methods into the shape of the wavefront for intermediate values of  and by characterizing this behavior, this model could be used to describe biological scenarios such as tumor growth, where  would represent the rate at which the tumor cells were able to degrade ECM in the surrounding environment.APPENDIX B: NUMERICAL METHODS Equations ( 12) and ( 13) are solved numerically subject to no flux boundary conditions ( 14) in  at  = 0 and  =  using the method of lines on the one-dimensional spatial domain [0, ] where  > 0 is chosen to be sufficiently large to remove boundary effects.In most cases, we take  = 200.The spatial domain is uniformly discretized with spacing Δ = 0.1 between each of the  = 1, … ,  spatial points, and the following discretization is used 43 : where   represents the value of  at the spatial point .For the model ( 12) and ( 13), we use this discretization twice, with  = (1 − ),  =  and for the second term in the flux as  = ,  = .Equations ( 12) and ( 13) can then be rewritten as a system of 2 ordinary differential equations given by for 1 ≤  ≤  − 1.To implement the boundary conditions, we introduce the ghost points  −1 and  +1 44 and set  0 () =  −1 (),  +1 () =   (), ∀ ≥ 0, (B4) so that We solve the system of equations (B2) and (B3) and (B5) and (B6) using the built-in Python solver scipy.integrate.solve_ivp with the explicit Runge-Kutta integration method of order 5 and time step  = 1.Convergence checks were completed by considering a range of tolerances, time, and spatial steps to ensure that the parameters used for simulations produced solutions within the second-order error associated with the numerical scheme.

APPENDIX C: COMPARISON TO OTHER MODELS IN THE LITERATURE
This study focuses on the impact of introducing volume-filling effects of cells and ECM to a model of cell invasion into ECM.There are a number of PDE model simplifications in the literature, including the following model, proposed as a minimal model for tumor growth into ECM in Colson et al. 15 : that assumes cell motility to be impacted by the presence of surrounding ECM only and cell proliferation impacted only by other cells, that is, the resource limiting cell proliferation is not space.Another similar model is presented in Browning et al. 35 to describe melanoma growth into skin and it is subsequently analyzed in El-Hachem et al. 14 The model can be interpreted to assume that cell motility is decreased by ECM and that cell proliferation is impacted by both other cells and ECM: The model variables and parameters are interpreted in the same way as in the model presented in this work (12) and (13).
We are particularly interested in comparing the population-level behaviors of the PDE model for cell invasion into ECM presented in this work, which incorporates volume-filling effects into both diffusion and proliferation of cells, to the simpler models without these volume-filling effects, presented in the literature.By looking at Figure C1, we can draw the following conclusions: all three models produce traveling wave solutions with a speed  ≥  min , where  min is the minimum speed predicted by standard traveling wave analysis.In fact, all of these speeds are dependent on both the initial density of ECM ahead of the wave,  0 , and the rescaled ECM degradation rate .The two models with the same reaction (growth) terms, depending on both cell and F I G U R E C 1 The relationship between the numerically estimated speed of traveling wave solutions to the system (C1) and (C2) on the left (blue), (C3) and (C4) in the middle (green), and ( 12) and ( 13) on the right (red), subject to the initial conditions ( 16) and (17).The numerically estimated traveling wave speed is obtained by tracing the point () such that ((), ) = 0.1.Further specifics of the parameter values and the numerical methods used can be found in Appendix B.
ECM preventing growth, predict the same traveling wave speed  min = 2(1 −  0 ), that is achieved numerically for  → 0 + .However, the model (C1) and (C2) presented in Colson et al. 15 predicts a speed  min = 2 √ 1 −  0 , that is also revealed for  → 0 + .As a result, the behaviors observed in these models for small rescaled ECM degradation rates  can be reproduced by studying a FKPP model (11) with the appropriate parameters.In the same manner, by looking at Figure C1, it is clear that all three models produce traveling waves with speed  → 2 − as  → ∞.The behaviors observed here can be studied by considering the standard FKPP model ( 18) with all parameters equal to unity.The FKPP model ( 18) is also a suitable model simplification for all three systems when  0 = 0.
The transition between the two asymptotic regions is yet to be fully characterized for any of the models, but it is clear that  is a monotonic, increasing function of  and  0 for all of the models.The critical value above which  begins to influence the speed is similar across the models but clearly depends on  0 and takes larger values across the models as more volume-filling effects are taken into account.Following intuition, we also find that, in general, the speed of invasion is slower as volume-filling effects are considered to impact more aspects of cell behaviors (from left to right in Figure C1).
The most obvious difference between these results is that the model ( 12) and ( 13) derived in this work does not permit traveling waves for compactly supported initial conditions in  when  0 = 1.This is a direct result of consistently including volume-filling effects across all the mechanisms of cell movement, such that there is always a maximum number of cells present at any point in space.The results match those of the model (C3) and (C4) when noncompactly supported initial conditions are simulated, as presented in Appendix A.
As such, the model presented in Browning et al. 35 provides a good model simplification by which to study the qualitative properties of the solutions to ( 12) and ( 13) across all parameter values when  0 ∈ [0, 1), with simplifications to the FKPP model also being appropriate as  → 0 + and  → ∞. so that the Jacobian is given by  12) and ( 13) subject to the initial conditions ( 16) and ( 17), for  0 = 0.2 in the top row and  0 = 0.8 in the bottom row and for rescaled ECM degradation rates  = 10 4 , 10 5 , 10 6 .Cell densities are shown in purple and ECM densities in orange.Further specifics of the parameter values and the numerical methods used can be found in Appendix B.
), for  0 = 0.2 in the top row and  0 = 0.8 in the bottom row, and for rescaled ECM degradation rates  = 5, 50, 500.Cell densities are shown in purple and ECM densities in orange at times  = 25, 50, 75, 100 from left to right.Further specifics of the parameter values and the numerical methods used can be found in Appendix B.
), for  0 = 0.2 in the top row and  0 = 0.8 in the bottom row, and for rescaled ECM degradation rates  = 10 −3 , 10 −2 , 10 −1 .Cell densities are shown in purple and ECM densities in orange at times  = 2500, 5000, 7500, 10000 from left to right.Further specifics of the parameter values and the numerical methods used can be found in Appendix B.
), for  0 = 0.2 in the top row and  0 = 0.8 in the bottom row, and for rescaled ECM degradation rates  = 10 4 , 10 5 , 10 6 .Cell densities are shown in purple and ECM densities in orange at times  = 25, 50, 75, 100 from left to right.Further specifics of the parameter values and the numerical methods used can be found in Appendix B.

A
C K N O W L E D G M E N T S R. M. C. is supported by funding from the Engineering and Physical Sciences Research Council (EPSRC) and Wolfson College, University of Oxford.T. L. gratefully acknowledges support from the Italian Ministry of University and Research (MUR) through the grant "Dipartimenti di Eccellenza 2018-2022" (Project no.E11G18000350001), the PRIN 2020 project (No. 2020JLWP23) "Integrated Mathematical Approaches to Socio-Epidemiological Dynamics" (CUP: E15F21005420006), and the INdAM group GNFM.The authors are grateful to Kevin Painter and Chloe Colson for interesting discussions regarding traveling waves in cell invasion models.D ATA AVA I L A B I L I T Y S TAT E M E N T Data sharing is not applicable to this article as no datasets were generated or analyzed during the current study.O R C I D Rebecca M. Crossley https://orcid.org/0000-0001-7342-0207F I G U R E A 1 Numerical solutions to the system (12) and (13) subject to the initial conditions (A1) and (A2) (panel (A)) or (A1) and (A3) (panel (B)), for  0 = 1 and  = 250.Cell densities are shown in purple and ECM densities in orange at times  = 2, 4, 6, 8, 10, 12, 14, 16 (from left to right) in panel (A) and times  = 25, 50, 75, 100 (from left to right) in panel (B).Note that the axis in the plot in panel (a) is zoomed in on  ∈ [0, 5] to display the initial behavior in the transient region before invasion stops.Further specifics of the parameter values and the numerical methods used for simulation can be found in Appendix B.