Travelling-Wave and Asymptotic Analysis of a Multiphase Moving Boundary Model for Engineered Tissue Growth

We derive a multiphase, moving boundary model to represent the development of tissue in vitro in a porous tissue engineering scaffold. We consider a cell, extra-cellular liquid and a rigid scaffold phase, and adopt Darcy’s law to relate the velocity of the cell and liquid phases to their respective pressures. Cell–cell and cell–scaffold interactions which can drive cellular motion are accounted for by utilising relevant constitutive assumptions for the pressure in the cell phase. We reduce the model to a nonlinear reaction–diffusion equation for the cell phase, coupled to a moving boundary condition for the tissue edge, the diffusivity being dependent on the cell and scaffold volume fractions, cell and liquid viscosities and parameters that relate to cellular motion. Numerical simulations reveal that the reduced model admits three regimes for the evolution of the tissue edge at large time: linear, logarithmic and stationary. Employing travelling-wave and asymptotic analysis, we characterise these regimes in terms of parameters related to cellular production and motion. The results of our investigation allow us to suggest optimal values for the governing parameters, so as to stimulate tissue growth in an engineering scaffold.


Introduction
In vitro tissue engineering is a form of regenerative medicine which often involves seeding cells into a porous bio-engineered scaffold that allows nutrient transport, structural support and a means for cell signalling activity (Chan and Leong 2008). Subject to the correct environment and growth factors, the cells will develop into a functional construct that can be used to restore damaged tissues and improve concepts in pharmaceutical research such as experimental drug therapy (Jensen et al. 2018 contributions from an assortment of scientific fields, tissue engineering is considered an interdisciplinary practice that has the potential to benefit a substantial proportion of the global population with devastating soft tissue, bone and whole organ diseases (Dzobo et al. 2018). The field of tissue engineering has enjoyed many successes; for example, the generation, replacement and longevity of engineered bones and bronchial tubes derived from the recipients' cells (Sato et al. 2008;Petite et al. 2000;Schimming and Schmelzeisen 2010). However, a shortage in the supply of donor tissue creates a demand on the field to make engineered tissue routinely clinically available (Levitt 2015). Whilst the field is rich in both theoretical and experimental knowledge, a lack of understanding regarding the processes by which cells assemble into tissues means that viable replacement constructs are only available in a minority of cases.
Extensive mathematical research has been undertaken to make sense of the complicated mechanisms within tissue growth. Some authors adopt a microscale approach, which can take the form of cellular automaton systems (Vitvitsky 2014;Lehotzky and Zupanc 2019;Youssef 2015) that seek to model interactions between a large number of individual cells. Whilst such systems can track the behaviour of cells, they can become computationally infeasible for tissue-scale simulations (Ermentrout and Edelstein-Keshet 1993). Some authors adopt a probabilistic approach (Fadai et al. 2019;Browning et al. 2019a, b) and exploit experimental data to explore ways in which tissue engineering techniques can be improved. For example, Sogutlu and Koc (2013) present a stochastic model to determine the expected required number of pores for each region of a tissue engineering scaffold.
Conversely, continuum models (O'Dea et al. 2010;Lemon et al. 2006;Breward et al. 2002;Eyles et al. 2019;Byrne et al. 2003 ) track the evolution of tissue constituents by employing systems of partial differential equations. Whilst continuum models cannot impose rules on individual cells, they can be derived by imposing averaging techniques (Drew 1983) on equations that govern cellular behaviour at a microscopic level, from which relevant mathematical techniques may be exploited to determine relationships between parameters and model behaviours. The study of in vitro tissue growth via continuum models is extensive, see O'Dea et al. (2012) and Klika et al. (2016) for reviews; of particular relevance to this study, however, Lemon et al. (2006) considers a continuum multiphase model to investigate how mechanical pressures within growing tissue influence the aggregation or dispersion of cells in a scaffold, and relates the existence of these regimes to the governing parameters. Lemon and King (2007a) examines travelling-wave solutions of the multiphase model formulated in Lemon et al. (2006) and find that in certain limits, the tissue propagates through the scaffold at a constant speed as either a forward or backward travelling wave, dependent on parameter values.
In this paper, we develop and analyse a continuum multiphase model that represents the development of tissue in vitro in an artificial scaffold. In our model, we aim to capture key features of tissue growth and extinction whilst developing a tractable formulation. In particular, we consider a porous flow description comprising a tissue cell phase, extra-cellular liquid phase and a scaffold phase, the former two being modelled as incompressible fluids and the latter as an inert solid. Tissue mechanics are accounted for by considering relevant constitutive assumptions in a similar fashion to those presented in Lemon et al. (2006) and Lemon and King (2007a). The model is reduced to a reaction-diffusion equation for the cell phase and a moving boundary condition for the tissue edge, after which travelling-wave, asymptotic and numerical methods are employed to deduce the resulting solution behaviour. The paper is constructed as follows. In Sect. 2, we formulate and subsequently reduce and non-dimensionalise the model. In Sect. 3, we present and discuss numerical solutions to the reduced model, which motivate the travelling-wave and asymptotic analyses conducted in Sects. 4 and 5. In Sect. 6, we draw some conclusions regarding the behaviour of the model and interpret the mathematical results in terms of the biological application.

Model Development
We construct a multiphase model to describe the growth of a nutrient rich tissue within a porous tissue engineering scaffold. For simplicity, we formulate the model in a onedimensional Cartesian geometry. The model consists of three phases: two of which are fluid phases denoted by n(x, t) and w(x, t), and represent the volume fraction of cells and extra-cellular liquid, respectively. A rigid, non-degradable scaffold with uniform volume fraction s is the third phase and remains constant, the porosity of the scaffold hence being given by 1 − s. Cell growth and death occur via mass transfer between n and w. The phases satisfy the no voids volume constraint: (1) The velocity fields v n (x, t) and v w (x, t), and pressures p n (x, t) and p w (x, t), are associated with the phases n and w accordingly. The spatial domain of the tissue evolves over time due to cellular motion, so we track it with a moving boundary, x = L(t). In the subsections that follow, we state equations that govern mass transfer between n and w, as well as provide constitutive assumptions for v n , v w , p n and p w suitable to describe tissue growth in a scaffold. We state necessary initial and boundary conditions for the variables and the moving boundary L(t), and simplify and non-dimensionalise the model.

Governing Equations
We assume that cells proliferate and assemble daughter cells from the available liquid and that when cells die, they decompose and dissolve into the liquid phase. In view of these processes, it is reasonable to follow Lemon et al. (2006), Byrne et al. (2003), Breward et al. (2002) and Preziosi and Tosin (2003) (and many others) and assume the densities of n and w to be equal. Following these assumptions, the mass transfer equations can be represented as where Γ is the net rate of cell proliferation. Adding the equations from (2) results in the overall conservation of mass condition where (1) has been used to eliminate the time derivative and to replace w with 1−n −s. Noting that n and w are modelled as fluids and s as a porous scaffold, we take the interphase drags to be dominated by those with the scaffold and neglect that between the tissue and liquid. In view of this, we apply Darcy's law to relate the velocity of the cell and liquid phases to their respective pressures. Following King and Franks (2004) and Eyles et al. (2019), we take where μ n and μ w represent the viscosity of the cell and liquid phases and K is the permeability of the scaffold. Remaining consistent with Lemon et al. (2006) and Lemon and King (2007a, b), we relate the cellular and extra-cellular liquid pressures via where represents extra pressures that arise due to cell-cell and cell-scaffold interactions. Since the scaffold is assumed to be inert and of uniform porosity, we suppress the dependence has on s from hereon for brevity. We note that combining (5) with the relations from (4) allows the elimination of p n and p w and provides

Initial and Boundary Conditions
Assuming the tissue to be symmetric about its centre ( Naturally, the cell volume fraction is identically zero at the edge of the tissue: The moving boundary L(t) moves with the cell velocity, hence The initial distribution of n and tissue boundary position, respectively, are denoted by n(x, 0) = n 0 (x) and

Model Reduction
We reduce the model to a reaction-diffusion equation and a moving boundary condition. Integrating (3) and applying the boundary conditions from (7) provides Here, we note μ n and μ w are assumed to be independent of n for simplicity. Substituting (11) into the first of (2) provides the reaction-diffusion equation: Combining (9) with (11) provides the moving boundary condition: where the boundary condition from (8) implies the boundary condition on v n from (7) becomes

Constitutive Assumptions
We now define constitutive assumptions for Γ and that are suitable to describe tissue growth in a rigid scaffold. We assume that daughter cells are constructed via mitosis using the available liquid and that when cells die via apoptosis, they dissolve into the liquid. Thus, we take where r m and r a are the positive constant rates of cell mitosis and apoptosis, and (1) is used to replace w with 1 − n − s. Following Lemon et al. (2006) and Lemon and King (2007a, b), an appropriate expression for (n) is for ν ∈ R and positive constants δ n , δ s and χ. The first term in (16) represents repulsive forces exerted between the cells at high volume fractions, as characterised by the singularity at n = 1 − s. The second term represents the propensity for cells to disperse or aggregate, with ν taking a positive or negative value accordingly. The third term represents repulsive forces that occur due to cell-scaffold interactions, whilst the fourth describes attractive forces between the cells and scaffold. For simplicity, we take δ := δ n = δ s . We note that Φ(n) must be strictly positive to prevent negative diffusion in (12) and nonlinear degeneracy in (13). This is achieved when ν > 0, which is henceforth assumed. Physically, this corresponds to a tendency for cells to spread through the scaffold (Lemon et al. 2006).
The parameter κ is shown in subsequent sections to be of crucial importance to the qualitative features of the model solutions. Physically, κ represents the difference between the scaffold porosity and the ratio between the cell death and growth rates. We note that the scaffold permeability parameter K , as seen in (11) is not present in the non-dimensional model (18)-(20). However, the scalings selected in (17) imply that the dimensional tissue boundary position increases with the scaffold permeability. A linear stability analysis around the steady states of (18), n = 0 and n = κ, is conducted in "Appendix A" and shows that the former is stable when κ < 0 and the latter when κ > 0. In view of this, we are primarily motivated to investigate (18)-(20) for different values of κ, though variations in s, μ, and η will also be considered in part so as to deduce their optimal values for the stimulation of tissue growth. Unless otherwise stated, we take μ = η = 1 and we adopt the initial conditions so that ω denotes the cell volume fraction at x = 0. Following Lemon and King (2007a, b), and unless otherwise stated, we set s = 0.2 and ω = 0.03, the former corresponding to a scaffold with a porosity of 0.8 and is consistent with the experimental study presented in Malda et al. (2004).

Numerical Results
We present and discuss the numerical solutions for n(x, t) and L(t) from the PDE system (18)-(20), paying separate attention to the cases κ > 0, κ < 0 and κ = 0. For numerical convenience, we fix the moving boundary by introducing the variable transform ξ = x/L(t) so that ξ ∈ [0, 1], which means (18) and (19) become Subject to the transformed boundary and initial conditions from (20), we numerically integrate (23) and (24) by discretising first-and second-order spatial derivatives using second-order finite differences. Upwind finite differences were used for the second term of (23). Temporal derivatives are numerically integrated by utilising ode23s in MATLAB.
For κ = 0.3, as shown in Fig. 1a, b, we observe semi-infinite travelling waves in n and linear growth in L after a period of transient growth from their initial states. For κ = 0, as shown in Fig. 1c, we observe n decaying from the initial data. Figure 1d shows unbounded growth in L. The inset shows L(t) and the function ln(t)/ √ 2 − 1 plotted against ln(t), from which we conclude that L grows logarithmically at large time. For κ = −0.3, as shown in Fig. 1e, we observe that n decays from the initial data more quickly than for κ = 0. The initial growth of L shown in Fig. 1f occurs due to the diffusion of n from the initial state; however, we observe the eventual formation of a steady state. Numerical simulations that are not included here suggest that travellingwave and steady-state behaviour is exhibited by (18)-(20) for all κ > 0 and κ < 0 accordingly.
Clearly, the case in which κ > 0 corresponds to effective tissue growth. This motivates a travelling-wave analysis of (18)-(20) for κ > 0 which is presented in Sect. 4, where we express the speed of the tissue edge in terms of the governing parameters. In Sect. 5, asymptotic solutions for n and L are found when 0 < κ 1, so that the cell distribution and tissue speed are explicitly available. Whilst the case κ < 0 results in tissue decay, an asymptotic analysis of (18)-(20) for this case is presented in Sect. 5. Overall, the results in this section suggest that κ > 0 must hold for tissue growth to occur, thus suggesting that tissue engineers should ensure that the porosity of the scaffold is at least larger than the ratio between the rate of cell death and growth. 4 Travelling-Wave Analysis for Ä > 0 Figure 1a, b indicates the emergence of semi-infinite travelling waves of constant speed for κ > 0. In light of this, we assume that for sufficiently large time, L ∼ ct where c is the constant wave speed at which the tissue edge moves. In this section, we employ travelling-wave analysis to obtain the wave speed c in terms of the governing parameters when κ > 0.
We transform (18)-(20) via the travelling-wave coordinates Following Fadai and Simpson (2020) and Fadai (2021) , we define Multiplying (25) by φ(n) and rewriting the conditions from (26) in terms of q(z), we obtain Here, we note the second boundary condition from (26) transforms into the second of (29) because φ(0) = 1. Dividing (28) by (27) we have Using the shooting method (presented in "Appendix B") to find the heteroclinic connection q(n) that connects (n, q) = (κ, 0) to (0, −c), we can determine a numerical approximation of the wave speed in terms of the governing parameters κ, s, μ and η. In Fig. 2a, b, the solid black line represents the relationships c(κ) and c(s), respectively, when (30) and (31)  line represents these wave speeds when obtained by numerically solving (18)-(20), and computing c by evaluating dL/dt at large time. In view of the close agreement between these two approaches to computing c, we henceforth concentrate on solutions provided by (30) and (31) for simplicity.
The results presented in Fig. 2 suggest that larger κ and smaller s increase the speed at which the tissue front grows, and since κ = 1 − s − r a /r m , this further corresponds to minimising r a and maximising r m and the porosity of the scaffold. In Fig. 3a, b we present the wave speeds c(μ, η) for κ = 0.3 and κ = 0.7, respectively. These results suggest that, for a fixed value of κ, the wave speed is maximised when μ, η → 0. Physically, this corresponds to the case where the viscosity of the cells is much greater than the viscosity of the liquid, and where repulsive forces exerted due to cell-cell and cell-scaffold interactions at high cell volume fractions dominate inter-cellular forces that give rise to cell dispersal. Furthermore, Fig. 3a, b indicates that the dependence c has on μ is weaker for κ = 0.3 than κ = 0.7. This suggests that for smaller κ, cell-cell and cell-scaffold interactions which can drive cellular motion are more prominent in controlling the wave speed than the cell and liquid viscosities. Additionally, and in agreement with Fig. 2a, Fig. 3 indicates that the wave speed increases as κ increases.

Asymptotic Analysis for |Ä| 1
In this section, we construct asymptotic solutions for n(x, t) and L(t) for |κ| 1 when t 1. Since scaffold porosity is a readily controllable parameter (in contrast to cell growth and death), the analysis in this section when κ = 1 − s − r a /r m > 0 can be associated with the case in which the scaffold porosity is low and tissue growth is successful, so that 0 < κ 1. The numerical results in Fig. 1e indicate that n 1 holds for κ < 0 when t 1. Furthermore, given that max{n} = κ when κ > 0, n 1 and t 1 is also expected when 0 < κ 1 for sufficiently large time. We therefore have φ(n) ∼ φ(0) = 1 when |κ| 1 and t 1. To aid the subsequent asymptotic analysis, we introduce which imply that (18) and (19) are simplified to To analyse the behaviour of (33) and (34), we follow Newman (1980) and adopt the ansatz wherein 0 < B < A 1. Imposing N (L, T ) = 0 on (35), we obtain that In view of (35), (33) provides In this section, initial conditions for N and L are chosen to satisfy (35) and (36) where α = 2B −2 0 (A 2 0 − B 2 0 ) 3/2 and β = sinh( We now deduce the large-T behaviour of N (x, T ) and L(T ), from which the largetime behaviour of n(x, t) and L(t) when |κ| 1 can subsequently be determined. Guided by the numerical results from Fig. 1d, the evolution of L satisfies L L 0 for sufficiently large T , so we have from (39) that which is then inverted to give for T 1. Equation (38) and the leading-order term in the above expansion are used to find the following large-T approximation for N : We now exploit (42) and (41) to deduce the large-time behaviour of n(x, t) and L(t) when |κ| 1.

Large-Time Behaviour of n and L When |Ä| 1
When 0 < κ 1, then n(x, t) takes the form of a travelling wave of constant speed and L L 0 for t 1. In contrast, when κ < 0, the numerical results in Fig. 1f suggest that L → L ∞ as t → ∞ for some finite constant L ∞ . In general, a large-time solution for L is unavailable given that L ∞ L 0 does not necessarily hold when t 1; however, when |κ| 1, then t = T to leading order at t = O(1), and L evolves according to (41) until t = O(1/|κ|). Since L L 0 when |κ| 1 and t 1, for |κ| 1 and t 1. Equation (43) implies that L ∼ κt/ √ 2 when 0 < κ 1 and t 1. Therefore, travelling waves propagate with speed c ∼ κ/ √ 2 when 0 < κ 1, this being in agreement with the numerical results from Fig. 2a. For 0 < κ 1, (43) implies that the growth of the tissue edge is logarithmic until t = O(1/κ) and linear thereafter. If tissue growth is successful, this suggests the formation of travelling waves with constant speed is delayed when the scaffold porosity is low. When κ < 0, the exponential terms from (43) are negligible as t → ∞ and we obtain when κ < 0 and |κ| 1. The leading order logarithmic term in (43) and (38) is used to find the following large-time approximations for n when |κ| 1 : We note that (43) and (45)  In Figs. 4 and 5, we compare the numerical solution for n and L when obtained by numerically solving the PDE system from (18)-(20) for κ = 0.001 and κ = −0.001 against their respective asymptotic solutions from (45), (46) and (43). Overall, an excellent agreement between the numerical and asymptotic solutions is observed. The large-T behaviors of N and L characterised by (41) and (42), and hence the asymptotic approximations from this subsection, are only valid for initial conditions that satisfy (35) and (36). We now show that solutions of (33) and (34) converge to solutions similar to that of (41) and (42) for a wider class of initial data.

Convergence of Asymptotic Solutions
Since the choice of initial cell distribution within the scaffold is likely to vary substantially in practice, it is important to determine the large-T behaviour of N and L for a wider class of initial data, such as those from (22). The asymptotic behaviour of (33) as T → ∞ comprises an interior layer near the interface within which the similarity reduction Following familiar arguments to that of the porous-Fisher equation (Aaronson 1980;Murray 2002), although we emphasise that N is not a travelling wave of the usual form, the solution to (47) is given by implies that the solution to which that matches into the exponential terms in (48), and hence the corresponding term in the interior layer, is given by for some unknown constant m, and this dominates the asymptotic behaviour of (50) as T → ∞. Therefore, (49) becomes where m depends on the initial data. Since N (L, t) = 0, (52) implies that for T 1. By comparing (52) and (42), we see that the asymptotic structure is retained despite the initial cell distribution for large T . In addition, (53) suggests that the choice of initial cell distribution does not affect the speed at which the tissue edge moves for large T , but does affect the position of the tissue boundary. We note that if N (x, 0) is chosen to satisfy (35), then (42) indicates that m = √ 2/α. In Fig. 6, we compare the numerical solutions for N and L when obtain by numerically solving the PDE system from (33) and (34) against the asymptotic solutions from (52) and (53) for two choices of N (x, 0) and L 0 = 1. The value of L 0 was found by solving N (L 0 , 0) = 0. For both N (x, 0), we are able to choose an m that provides excellent agreement between the numerical and asymptotic solutions. We note that the large-time behaviour for n(x, t) and L(t) when κ = 0 can be extracted directly from (52) and (53) given that N = n when κ = 0 and lim κ→0 T = t. This justifies the numerical results observed in Fig. 1c, d. extra-cellular liquid phase and a scaffold phase and adopt Darcy's law to relate the velocity of the cell and liquid phases to their respective pressures. The model includes mechanisms to represent cell growth and death, and pressures that arise from cell-cell and cell-scaffold interactions. We employ a moving boundary, x = L(t), to track the speed at which the tissue edge propagates through the scaffold. We reduce the model to a nonlinear reaction-diffusion equation for the cell volume fraction, n(x, t), and a moving boundary condition for the tissue edge. The diffusivity of the reactiondiffusion equation is dependent on the cell and scaffold volume fractions; cell and liquid viscosities, and pressures that arise from cell-cell and cell-scaffold interactions. Non-dimensionalisation of the model shows that the tissue boundary position increases with the scaffold permeability and exposes important dimensionless groupings. One such grouping, κ, that describes the difference between the scaffold porosity and the ratio between the cell death and growth rates is of crucial importance to the qualitative features of the cell phase evolution. The model admits three regimes for the evolution of the cell volume fraction and the moving boundary, based on the sign of κ. Employing travelling-wave and asymptotic analysis, we characterise these regimes in terms of κ and parameters related to cellular motion.
The case in which κ > 0 corresponds to the successful growth of tissue, which suggests that tissue engineers should ensure that the porosity of the scaffold is at least larger than the ratio between the rate of cell death and growth. For κ > 0, we show that the cell volume fraction, n(x, t), spreads through the scaffold as a semi-infinite travelling wave with constant speed, emerging from the steady state n = κ. Employing travelling-wave analysis, we accurately compute the wave speed (i.e., the speed at which the tissue edge moves through the scaffold) as a function of the governing parameters. We find that the wave speed is greatest when the rate of apoptosis is negligible in comparison with that of mitosis, the viscosity of the cells is much greater than the viscosity of the liquid, and when repulsive forces exerted due to cell-cell and cell-scaffold interactions at high cell volume fractions dominate inter-cellular forces that give rise to cell dispersal. We also find that the wave speed increases as the scaffold porosity increases; however, we note that the cells will require a sufficient amount of scaffold on which to attach, so an upper bound on the porosity is to be expected. Furthermore, we deduce that for smaller values of κ, and hence scaffolds with small porosities, cell-cell and cell-scaffold interactions which can drive cellular motion are more prominent in controlling the wave speed than the cell and liquid viscosities.
For |κ| 1, we employ asymptotic analysis to find explicit solutions for n and L. Since scaffold porosity is a readily controllable parameter (in contrast to cell growth and death), the analysis in this section when κ = 1 − s − r a /r m > 0 can be associated with the case in which the scaffold porosity is low and tissue growth is successful, so that 0 < κ 1. When |κ| 1, the growth of the tissue edge is logarithmic until t = O(1/κ) and linear thereafter, thus suggesting the formation of travelling waves with constant speed is delayed as κ → 0 + , and hence when the scaffold porosity is low. For κ < 0, we deduce that the cell volume fraction decays exponentially with rate κ at large time, with the moving boundary tending towards a steady state. For κ < 0 and |κ| 1, the evolution of the L is shown to be logarithmic until t = O(1/κ) and approaches a steady state thereafter, the value of which is found explicitly and related to κ and the initial conditions employed in the model. For |κ| 1, we also demonstrated that the choice of initial cell distribution does not affect the eventual distribution of cells within the scaffold, nor the speed at which the tissue edge moves, but does affect the position of the tissue boundary. For a functional tissue construct to develop within a scaffold, cells must be exposed to the correct environment and stimulated with growth factors such as oxygen. There must also be a sufficient amount of scaffold on which the cells can adhere. Whilst key features of tissue growth such as cell mitosis, apoptosis and motion are included in this paper, concepts such as environmental pressures, cellular adhesion, and nutrient supply have not been considered. Therefore, following Lemon and King (2007a), a natural extension of this work would include examining the influence that nutrient limitation has on cell growth. We leave these extensions for future consideration.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

A Linear Stability Analysis
A linear stability analysis around the steady states of (18), n ∞ = 0, κ, provides insight into the dependence of the model behaviour on κ and φ(n). Neglecting the influence of the moving boundary condition on the stability of (18), we linearise on a semi-infinite domain. We substitute n = n ∞ + ε exp(iγ x + λt) into (18) for a perturbation wave number γ and growth rate λ where ε 1. Considering terms of O(ε) only, the growth rate for perturbations of wave length 2π/γ is λ = −n ∞ γ 2 φ(n ∞ ) + 2 + κ.
Since φ(n) is assumed to be positive for any n, the steady state n ∞ = κ is stable for all κ > 0. For n ∞ = 0, we have λ = κ which indicates stability when κ < 0.

B The Shooting Method
In this section, we formulate a numerical shooting method to find the wave speed c, as stated in Sect. 4. To do this, we send trajectories of q(n) from (κ, 0) to find a trajectory that connects to (0, −c). Computationally, it is more straightforward to shoot trajectories forwards as opposed to backwards, and we hence introduce the change of variable X = κ − n and obtain Noting that the end points are computationally singular, we use q(ζ ) = −ζ and q(κ − ζ ) = −c where 0 < ζ 1 is some user-defined tolerance. We employ the discrepancy function where q c (κ − ζ ) is the solution to (55) evaluated at X = κ − ζ for a trial wave speed c. Equation (55) is numerically integrated with ode23s in MATLAB with the initial condition q(ζ ) = −ζ. Using fzero in MATLAB to find the zero of E(c), the wave speed c and the corresponding heteroclinic trajectory q(X ) is determined.