Bifurcations and dynamics emergent from lattice and continuum models of bioactive porous media

TL;DR

This paper compares lattice and continuum models of bioactive porous media, revealing that lattice models exhibit asymmetric behaviors and bifurcations absent in continuum models, due to their structural differences.

Contribution

It demonstrates how the discrete lattice structure allows for complex asymmetric dynamics and bifurcations not present in the continuum model, advancing understanding of such systems.

Findings

Lattice models show asymmetric oscillations and steady states.

Continuum models exhibit only symmetric solutions.

Nonlocal reaction-diffusion mechanisms drive oscillations.

Abstract

We study dynamics emergent from a two-dimensional reaction--diffusion process modelled via a finite lattice dynamical system, as well as an analogous PDE system, involving spatially nonlocal interactions. These models govern the evolution of cells in a bioactive porous medium, with evolution of the local cell density depending on a coupled quasi--static fluid flow problem. We demonstrate differences emergent from the choice of a discrete lattice or a continuum for the spatial domain of such a process. We find long--time oscillations and steady states in cell density in both lattice and continuum models, but that the continuum model only exhibits solutions with vertical symmetry, independent of initial data, whereas the finite lattice admits asymmetric oscillations and steady states arising from symmetry-breaking bifurcations. We conjecture that it is the structure of the finite lattice which allows for more complicated asymmetric dynamics. Our analysis suggests that the origin of both types of oscillations is a nonlocal reaction-diffusion mechanism mediated by quasi-static fluid flow.

Equations

$$\sum _{j=1}^{n^2}{A}_{ij}{R}_{ij}^4(\widetilde{p_i}-\widetilde{p_j})=\{\begin{array}{ll}1-\widetilde{p_i},& 1\le i\le n,\\ 0,& n<i\le n^2-n,\\ -\widetilde{p_i},& n^2-n<i\le n^2,\end{array}$$

$$\frac{d{N}_i}{dt}=F_1(p_i)N_i(1-N_i)-F_2(p_i)N_i+\delta n^2\sum _{j=1}^{n^2}{A}_{ij}(N_j-N_i),1\le i\le n^2,$$

$$R_{ij}\equiv R(N_i,N_j)=1-\frac{\rho }{2}(N_i+N_j),$$

$$p_i=\widetilde{p_i}\frac{1}{Q},$$

$$Q=\sum _{i=1}^{n}1-{\tilde{p}}_i,$$

$$F_1(p_i)=1-\left(\frac{1}{2}\right)(\tanh [g(p_i-p_l)]+1),$$

$$F_2(p_i)=\left(\frac{1}{2}\right)(\tanh [g(p_i-p_l)]+1),$$

$$N_i(0)=N_{i0}.$$

$$\mathbf{u}=-k(N)\nabla \tilde{p},$$

$$\nabla \cdot \mathbf{u}=0,$$

$$\frac{\partial N}{\partial t}=F_1(p)N(1-N)-F_2(p)N+\delta {\nabla }^{2}N,$$

$$Q_c={\int }_0^1-k(N(0,\hat{y},t))\frac{\partial p}{\partial x}(0,\hat{y})d\hat{y},$$

$$p=\frac{\tilde{p}}{Q_c},$$

$$k(N)=(1-\rho N{)}^3,$$

$$F_1(p)=1-\left(\frac{1}{2}\right)(\tanh [g(p-p_c)]+1),$$

$$F_2(p)=\left(\frac{1}{2}\right)(\tanh [g(p-p_c)]+1),$$

$$\mathbf{u}\cdot \mathbf{n}=0\text{at }y=0,1,0\le x\le 1,$$

$$\tilde{p}=1\text{at }x=0,0\le y\le 1,\tilde{p}=0\text{at }x=1,0\le y\le 1,$$

$$\mathbf{n}\cdot \nabla N=0\text{for }\mathbf{x}\in \partial [0,1{]}^2,N(x,y,0)=N_0(x,y).$$

$$\sum _{j=1}^{n^2}{A}_{ij}{R}_{ij}^4(p_i-p_j)=\{\begin{array}{ll}\frac{1}{Q}-p_i,& 1\le i\le n,\\ 0,& n<i\le n^2-n,\\ -p_i,& n^2-n<i\le n^2,\end{array}$$

$$Q=\frac{n}{1+{\sum }_{j=1}^n{p}_j}.$$

$$\sum _{j=1}^{n^2}{A}_{ij}{R}_{ij}^4(p_i-p_j)=\{\begin{array}{ll}\frac{1}{n}\left(1+{\sum }_{j=1}^n{p}_j\right)-p_i,& 1\le i\le n,\\ 0,& n<i\le n^2-n,\\ -p_i,& n^2-n<i\le n^2.\end{array}$$

$$\sum _{j=1}^n{\hat{A}}_{ij}{R}^4({\hat{N}}_i,{\hat{N}}_j)({\hat{p}}_i-{\hat{p}}_j)=\{\begin{array}{ll}\frac{1}{n},& i=1,\\ 0,& 1<i<n,\\ -{\hat{p}}_n,& i=n,\end{array}$$

$${\hat{p}}_i=\frac{1}{n}\left(1+\sum _{k=i}^{n-1}{R}^{-4}({\hat{N}}_k,{\hat{N}}_{k+1})\right).$$

$$\frac{d\widehat{N_i}}{dt}=F_1({\hat{p}}_i)\widehat{N_i}(1-\widehat{N_i})-F_2({\hat{p}}_i)\widehat{N_i}+\delta n^2\sum _{j=1}^n{\hat{A}}_{ij}({\hat{N}}_j-{\hat{N}}_i),1\le i\le n,$$

$$\frac{\partial }{{\partial }_x}\left(k(x)\frac{\partial \tilde{p}}{{\partial }_x}\right)=0,$$

$$Q_c=-k(N(0))\frac{\partial \tilde{p}}{\partial x}(0),$$

$$p=\frac{\tilde{p}}{Q_c}.$$

$$\tilde{p}(x)=\frac{{\int }_x^1\frac{1}{k(N(\alpha ))}d\alpha }{{\int }_0^1\frac{1}{k(N(\alpha ))}d\alpha }.$$

$$p(x)={\int }_x^1\frac{1}{k(N(\alpha ))}d\alpha .$$

Full text

{history}

1 Introduction

We study a reaction–diffusion process on a two–dimensional domain modelled by a finite lattice dynamical system, which includes nonlocal coupling between nodes. This model arises as a mathematical caricature of a model developed in [Krause et al., 2017] representing fluid and cell interactions in a bioactive porous medium. While lattice–dynamical systems have been studied in a variety of different settings [Buzano & Golubitsky, 1983; Winterbottom, 2004; Wang et al., 2006; Kamei, 2009; Sattinger, 1980], we note that processes defined on finite lattices are less well–studied, and have been shown to exhibit unique dynamics due to the influence of boundaries [Chow et al., 1996; Gillis & Golubitsky, 1997; Golubitsky et al., 2004]. Additionally, our model also exhibits nonlocal coupling across the lattice, which is also of interest from the dynamical systems perspective [Gourley et al., 2001; Billingham, 2004; Hamel & Ryzhik, 2014]. We also study an analogous spatially continuous system, with the chief goal of comparing the dynamics between these two paradigms, and hence elucidating the roles of the finite lattice and the nonlocal coupling in our reaction–diffusion system.

The spatial domain in our models represents a porous tissue scaffold seeded with cells and perfused with fluid in order to stimulate cell growth and the formation of viable tissue, which is a primary goal of the field of tissue engineering. A major challenge in tissue engineering is to understand and exploit processes that occur on very different spatial and temporal scales in tissues and organs Vafai [2010]. Such scales permit the development of mathematical models which are able to help elucidate complex physical and biological processes relevant for tissue engineering O’Dea et al. [2012]; Van Blitterswijk & Thomsen [2008]. In addition to providing conceptual understanding, quantitative mathematical models can be used to predict and optimize experimental operating regimes which saves on costly and time-consuming experimental trials Geris [2013]. While continuum models for fluid flow in porous media can be justified for pore networks with many pores, corresponding either to large porous materials (such as oil reservoirs), or to materials with densely connected pore structures (such as sponges) Bear [1972], many tissue scaffolds used experimentally have pore sizes and scaffold geometries such that the number of pores is $O(10^{3})$ or less [Cox et al., 2015; German & Madihally, 2016; Loh & Choong, 2013; Vafai, 2010], suggesting that lattice models may be more appropriate than continuum models, at least in some settings.

While lattice approaches exist in related literature, such as spatial models of cell monolayers [Byrne & Drasdo, 2009], tumour vascularization [Anderson & Chaplain, 1998], biofilms [Thullner & Baveye, 2008], angiogenesis [McDougall, 2002; Scianna et al., 2013], and even tissue engineering [Barbotteau et al., 2003], there is still much work needed to understand the impact of finite pore networks on porous media undergoing structural changes due to cell growth; see McDougall & Sorbie [1997]; Blunt [2001]; O’Dea et al. [2012] for a more thorough comparison between discrete and continuous models. Krause et al. [2017] introduced two models of a two-dimensional porous tissue scaffold perfused with fluid. These models incorporated interactions between a viscous fluid phase representing a culture medium, and a cell phase modelling the biomass of proliferating cells and their extracellular matrix. Cells were allowed to diffuse locally throughout the domain, but would die if the local fluid shear stress became too large. Cell proliferation was assumed to block pores, reducing the local permeability of the scaffold, and hence introducing a nonlocal interaction between cells at different parts of the domain. Krause et al. [2017] modelled the system using a continuum approach involving a system of PDEs, and a lattice approach involving a spatially embedded network of ODEs, and discussed quantitative and qualitative differences between the modelling predictions arising from these approaches. Oscillatory solutions were only present for the lattice model, and were more pronounced for smaller pore networks.

Spatially discrete and continuous dynamical systems have been explored in a variety of related fields to model cells and tissues [Fraternali et al., 2014; Murisic et al., 2015; Matsiaka et al., 2018], providing further motivation to compare dynamics emergent from each modelling approach. Symmetry breaking bifurcations on lattices and other network structures have found relevance in fluid mechanics [Crawford & Knobloch, 1991], optics [Kevrekidis et al., 2005], mathematical physics [Bendix et al., 2009; Brazhnyi & Malomed, 2011], chaotic systems [Pikovsky & Grassberger, 1991; Rothos et al., 2002], and pattern formation Wolfrum [2012], to name a few examples. For a review of symmetry breaking dynamics on lattices, see Sattinger [1980]. Bifurcations due to lattice structures have been classified on a number of crystal lattices [Janssen & Tjon, 1983; Silber et al., 1992], and such results also extend to a number of other systems on lattices Buzano & Golubitsky [1983]; Winterbottom [2004]; Wang et al. [2006]; Kamei [2009], with Hopf bifurcations in particular appearing quite commonly Silber & Knobloch [1991]. We remark that many of these systems are formulated on infinite lattices, yet for practical purposes, models relevant to tissue engineering will likely be posed on finite lattices. This finiteness of the lattice results in the possibility of additional broken symmetries due to boundary effects Chow et al. [1996]; Gillis & Golubitsky [1997]; Golubitsky et al. [2004]. Therefore, depending on the specific functional structure of the dynamical system, and of the underlying spatial lattice, symmetry preserving or symmetry breaking bifurcations may be ubiquitous to the model dynamics. Additionally, Duncan et al. [2015] has demonstrated multistability in discrete stochastic chemical systems, whereas the continuum limit of these systems admits only one stable equilibrium state. We will demonstrate symmetry–breaking bifurcations and multistability in our lattice system which is not present in the spatially continuous analogue.

Motivated by the aforementioned results, in the present paper we systematically compare the dynamics emergent from lattice and continuum models of bioactive porous media. We modify the models presented in Krause et al. [2017] by assuming that cells are sensitive to high pressures, rather than fluid shear stress. This is motivated by computational modelling of cell death due to hydrostatic pressure Byrne & Drasdo [2009]; Byrne & Preziosi [2003]; Nessler et al. [2016], and by experimental work demonstrating cell death due to pressure Mañas & Mackey [2004]. These modified “fluid pressure” dependent models are more amenable to mathematical and numerical investigation, due to the simpler coupling between the fluid flow and growth processes compared with the fluid shear stress coupling used in [Krause et al., 2017] (which implicitly involved gradients in pressure, rather than values of the pressure itself).

The paper is organised as follows. In Section 2, we present our lattice and continuum models, and briefly discuss the underlying modelling assumptions and components. In Section 3, we present numerical solutions of the lattice and continuum models, demonstrating typical solution behaviours. In Section 4, we derive asymptotic solutions for the cell density appropriate for large values of the diffusion parameter, and show that these always tend to a steady spatially uniform state. We use these solutions to bound the region in parameter space where we anticipate non-equilibrium dynamics (e.g. oscillations). In Section 5, we describe the more complicated dynamics outside of this asymptotic parameter regime. We do this by presenting numerical bifurcation diagrams over two parameters for large and small lattices alongside the continuum model. Using these diagrams we broadly characterize the steady states and long-time oscillatory states by their symmetry properties. We present continuation results showing the existence of Hopf and pitchfork bifurcations that break the vertical symmetry of steady state solutions. Additionally, we reveal that a different kind of oscillation exists, which preserves vertical symmetry and is not due to a local bifurcation from a steady state. Finally, we give an overview of these results in Section 6, and discuss implications of our results.

2 Lattice and continuum models of a bioactive porous media

We consider a square lattice of $n$ nodes per side, and at each node define a pressure $p_{i}$ and a cell density $N_{i}$ for $1\leq i\leq n^{2}$. See Figure 1 for a representation of this lattice. We label nodal variables with single indices, counting from the bottom-left upwards so that $p_{1}$ is the pressure at the bottom-left node, and $p_{n^{2}}$ is the pressure at the top-right node. We define the adjacency matrix $A$ by $A_{ij}=A_{ji}=1$ if nodes $i$ and $j$ are connected, and $A_{ij}=0$ otherwise. We write the graph Laplacian as the matrix $L$ such that $\sum_{j}A_{ij}(v_{i}-v_{j})=\sum_{j}L_{ij}v_{j}$, and note that this formulation of $L$ accounts for no–flux boundary conditions at the edges of the lattice. See Newman [2010] (and references therein) for a discussion of networks in general and properties of the matrices $A$ and $L$.

We assume the fluid is viscous, Newtonian, and incompressible. We model the flow rate between nodes $i$ and $j$ via Poiseuille flow in a pipe of radius $R_{ij}$, and assume that $R_{ij}$ depends on the cell density at nodes $i$ and $j$. We assume that a constant fluid flow enters the domain from the left and leaves the domain along the right side of the lattice. We assume impermeable horizontal boundaries. We prescribe a unit pressure drop across the scaffold which will have an associated fluid flow rate $Q$ into the domain, and we use a pressure rescaling procedure described in Krause et al. [2017], based on work in Shakeel et al. [2013], to ensure that a unit flow rate enters the domain (this is the rescaled physical boundary condition, corresponding to a constant fluid flow into the domain). We assume that there is a critical pressure $p_{l}$ so that cells at node $i$ grow logistically for values of the pressure $p_{i}<p_{l}$, and die for values of $p_{i}>p_{l}$. We also assume that cells can diffuse between nodes. Our system of non-dimensional equations is then

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where the parameter $\delta$ represents cell diffusion, $\rho$ determines the effect of cell density on the pipe radius and hence the local fluid flow, and $g$ is a smoothing parameter. We also prescribe the initial data,

[TABLE]

Equations (1) represent conservation of mass enforced at each node, with a unit pressure drop imposed across the square domain. Equations (2) describe logistic cell growth and death modulated by the functions $F_{1}$ and $F_{2}$ of the pressure at that node, along with cell diffusion between nodes. Equation (3) describes how the presence of cells influences the pipe radii, and Equations (4)-(5) rescale the pressures to match the total fluid flow into the domain. Equations (6)-(7) describe functions of the pressure chosen to model a ‘proliferation’ and a ‘death’ response of the cells to low and high values of fluid pressure respectively. The use of hyperbolic tangent functions is a modelling choice which follows Shakeel et al. [2013], although similar choices such as step functions and mollified step functions have been used elsewhere Coletti et al. [2006].

The analogous dimensionless continuum model is

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $\delta$, $\rho$, $g$ are positive constants with the same meaning as before, $p_{c}$ is analogous to $p_{l}$, and $Q_{c}$ is analogous to $Q$. We prescribe the boundary and initial conditions

[TABLE]

[TABLE]

[TABLE]

Both models have no-flux (Neumann) conditions on all boundaries for the cell density, and on the horizontal boundaries for the fluid. To describe solutions with certain symmetries, we also consider spatially 1-D models without horizontal boundaries, as well as periodic conditions on both cell density and fluid flux along the horizontal boundaries. For the lattice, periodic boundary conditions are encoded in the adjacency matrix $A_{ij}$ by simply adding an edge between each node at the top and the bottom of the lattice.

We now consider simplifications of the models available when there is a vertical symmetry - that is, when the cell density distribution (and hence the pressure) is vertically symmetric. This is motivated by the presence of this symmetry in some numerical solutions from Section 3, and will provide a basis for understanding some mathematical aspects of these models.

2.1 1-D lattice model

We first simplify the lattice equations by substituting Equation (4) into Equation (1) to give

[TABLE]

where $Q$ is given in terms of $p_{i}$ from Equation (5) by

[TABLE]

Using this expression for $Q$, our system of fluid equations (20) becomes

[TABLE]

We solve the pressure equations when the solutions are vertically symmetric. We label the value at a representative node by $p_{i}\equiv\hat{p}_{k}$ where $k=\left\lceil\frac{i-1}{n}\right\rceil$, and $\left\lceil a\right\rceil$ denotes the smallest integer greater than $a$. So the first $n$ variables have the value $\hat{p}_{1}$, the next $n$ have the value $\hat{p}_{2}$, etc. We analogously denote the cell densities as $N_{i}=\hat{N}_{k}$. We rewrite Equation (21) using these variables as

[TABLE]

where we have used the simplification $(1+\sum_{j=1}^{n}p_{j})/n-p_{i}=(1+n\hat{p}_{1})/n-\hat{p}_{1}=1/n$. The adjacency matrix $\hat{A}_{ij}$ represents the path graph on $n$ vertices, so that $\hat{A}_{ij}=1$ if $i=j\pm 1$ or $\hat{A}_{ij}=0$ otherwise. Summing all $n$ equations in (22) we find that $\hat{p}_{n}=1/n$, where we have used the symmetry of the function $R$. For $1\leq i<n$ we can iteratively solve these equations to find

[TABLE]

This is equivalent to saying that resistance in serial circuits is additive, where this pressure variable is playing the role of voltage and $R^{-4}(\hat{N}_{k},\hat{N}_{k+1})$ is the resistance between nodes $k$ and $k+1$ [Oh et al., 2012].

Hence, in the 1-D model, we obtain a system of $n$ coupled equations for the cell density:

[TABLE]

where $\hat{p}_{i}$ given by (23).

2.2 1-D continuum model

We assume that the cell density varies only in the $x$ direction, so that the permeability $k=k(N(x))$ also only varies in this direction. We then have by Equations (9)-(10) and (17)-(18) that $\tilde{p}=\tilde{p}(x)$, so the governing flow equations become

[TABLE]

[TABLE]

[TABLE]

Integrating Equation (25) subject to the boundary conditions (18) we find

[TABLE]

Substituting (28) into (26) we compute $Q_{c}$, and from (27) we find

[TABLE]

Substituting (29) into Equation (11) we obtain a single equation for the cell density,

[TABLE]

Equation (30) is a spatially nonlocal scalar reaction-diffusion equation for the cell density $N$.

2.3 Boundedness of dynamics

We now comment on the boundedness of solutions to both models in the full 2-D setting, which will enable us to deduce reasonable ranges of model parameters n Sections 3 and 5. We first show that solutions to the Equations (1)-(7) are bounded within the region $[0,1]^{n^{2}}$. We note that for all $p\in\mathbb{R}$, the functions in Equations (6)-(7) satisfy $0\leq F_{1}(p)\leq 1$ and $0\leq F_{2}(p)\leq 1$. Now we assume that for some $t^{*}$, $N_{i}(t^{*})\in[0,1]$ for all $i$, and $N_{k}(t^{*})=1$ for at least one $k$. We then have that

[TABLE]

Hence continuously differentiable ($C^{1}$) solutions to Equations (2) are bounded above.

Similarly, assume that for some $t^{*}$, $N_{i}(t^{*})\in[0,1]$ for all $i$, and that $N_{m}(t^{*})=0$ for at least one $m$. We then have that,

[TABLE]

and so continuously differentiable solutions are also bounded below. Note that these bounds are independent of the pressures $p_{i}$. So we have that $N_{i}(t)\in[0,1]$ for all $i$ and for all $t\geq 0$.

While we are not aware of a formal comparison principle applicable to Equations (9)-(19), the above argument gives a heuristic bound on the behaviour of solutions for these equations as well. The existence of spatiotemporal oscillations, which we show in Section 3, suggests that a general comparison principle does not hold for this system. Hence a rigorous treatment of Equations (9)-(19) would be needed in order to formally extend this argument. Nevertheless, we will assume that PDE solutions are similarly bounded, and our numerical solutions confirm this.

3 Exploration of solution behaviours

In this section we present numerical solutions for the 2-D lattice equations (1)-(7) as well as the PDE system given by (9)-(19). We fix $\rho=0.9$ so that cell growth significantly affects the effective permeability of the medium. To capture sufficiently sharp behaviour in the pressure functions, we set $g=60$. The remaining parameters are varied to demonstrate the range of solution behaviours. The PDE system was simulated using the finite element software COMSOL with 24,912 triangular elements, as well as a finite difference scheme to verify the results. The lattice model was solved using an explicit adaptive Runge-Kutta method, ‘ode45’ in Matlab, with a maximal time step of $10^{-3}$. The restriction of the maximal time step was to ensure a good approximation of bifurcation phenomena Christodoulou [2008]. Convergence checks were carried out in the time step in the numerical schemes for both models, and in the number of elements for the PDE.

From Section 2.3 and the monotonicity of Equations (23) and (29), we can compute the maximum and minimum possible pressures throughout the domain in each (1-D) model. We then deduce the ranges of the pressure threshold parameters $p_{l}$ and $p_{c}$ that lie between maximum and minimum pressures in the scaffold. To compute these extremal pressures, we set the cell densities $N_{i}=0$ or $N_{i}=1$ for all $i$. Using (23), we see that if $p_{l}<1/n$, then $p_{i}>p_{l}$ for all $i$, so that cells at every node will die (exponentially decay) due to high values of the pressure. Similarly, if $p_{l}>(1+(n-1)(1-\rho)^{-4})/n$, then $p_{i}<p_{l}$, and all nodes will approach a uniform value of $N_{i}=1$ logistically. In order to realize nontrivial dynamics, we vary $p_{l}$ as

[TABLE]

for the lattice and using (29), similar reasoning leads to the bounds,

[TABLE]

for the PDE. Note the difference in the orders of magnitude between these two ranges, which is due to the constitutive difference between the response of the flow to the cell density. Specifically, the Poiseuille flow in the lattice given by (1) has a quartic dependence on cell density, whereas the cubic permeability in Darcy’s law given by (14) is cubic in the cell density. See Krause et al. [2017] for a discussion of why these constitutive exponents differ between the models. Changing the exponent in the permeability, Equation (14), from $3$ to $4$ gives quantitatively similar behaviour for the PDE and very large lattices.

For all of our simulations, we set the initial data to be a uniform state perturbed by normally distributed spatial noise with zero mean and standard deviation $10^{-4}$, so that $N_{i}(0)\approx 0.1$ and $N_{0}(x,y)\approx 0.1$ . The lattice and PDE models are both insensitive to perturbations in initial conditions, giving identical long-time behaviour across hundreds of realizations of these noisy initial data. The cell densities grow uniformly until the pressure at some point of the domain reaches the threshold $p_{c}$ for the PDE or $p_{l}$ for the lattice, where some regions then experience cell death due to high pressure. The simulations then tend to a spatial horizontal gradient in cell density with low cell density on the left side of the lattice (where the pressure is highest) and high cell density on the right side of the lattice (where the pressure is lowest). We observe one of three different long time behaviours: a stable steady state, a vertically symmetric ‘pulsating’ oscillation, or a vertically asymmetric oscillation. The steady states we observe are either vertically symmetric or asymmetric.

At very large values of $\delta$ (e.g. $\sim O(1)$), or large values of the threshold parameters $p_{l}$ or $p_{c}$, all simulations tend to a similar steady state with a gradual horizontal gradient in cell density (not shown). In Figure 2 we plot two example cell density distributions alongside corresponding plots of pressure for $\delta=10^{-1}$. We also plot time series of the nodal evolution along the diagonals of the scaffold, i.e. all nodes along the lines $x=y$ and $x=1-y$, in Figures 2e-2f. For large $\delta$, there is a non-negligible cell density present in the region of high pressure despite the local cell density being governed by exponential death. The pressure does not vary significantly in this region of lower cell density near the left side of the scaffold, as the effective permeability changes little at low cell densities. In these parameter ranges, the 1-D reduction described in Section 2.1 is a quantitatively accurate approximation to the cell density along any row of the lattice.

For smaller values of the diffusion we observe a different long-time behaviour. In some cases, the cell density throughout the scaffold oscillates with a ‘pulsing’ between the right and the left side of the domain, maintaining a vertical symmetry throughout the oscillations. In Figure 3 we show the cell densities across diagonal nodes for two cases: a lattice of size $n=100$ in Figure 3a and the PDE in Figure 3b. For each, we observe stable oscillations throughout the domain. The lattice simulation in Figure 3a has the largest oscillations in cell density for nodes at the far left of the lattice, and these values never get close to the zero. In comparison, the PDE simulation in Figure 3b has its largest oscillations in the middle of the domain, and the far leftmost part of the domain exhibits very small cell densities. In Figure 4 we plot the cell density at three different points in time during one oscillation half period for the same simulations as in 3. In Figure 4a, we plot just the 1-D values of the column-averaged cell density as these oscillations are vertically symmetric, and similarly we plot vertically-averaged cell densities for the PDE in Figure 4b. The oscillations move between the yellow and blue curves with a period of $t\approx 0.94$ in Figure 4a and $t\approx 1.12$ in Figure 4b. These large oscillations persist in the 1-D reductions of both models, as well as when the horizontal boundary conditions for cells and fluid flow are changed from Neumann to periodic conditions. This suggests that these oscillations are inherently a 1-D phenomenon.

The aforementioned behaviours - vertically symmetric steady states and oscillations - appear in the PDE and in large lattices (e.g. for $n=100$). Qualitatively different steady states and oscillations exist in the smaller lattices. In many of the smaller lattice simulations, we observe a complicated breaking of the vertical symmetry. We give an example of vertically asymmetric oscillations in Figure 5, where the cell density oscillates vertically between the two impermeable horizontal boundaries. We plot a corresponding time series in Figure 6a, which demonstrates a separation of solutions into ‘bands,’ corresponding to specific columns in the lattice. The largest of these bands that changes between a growing and a dying state (shown in dashed red and solid green in Figure 6a) corresponds to the third column from the right boundary in the lattice. The bands which are always green (locally growing logistically) correspond to nodes to the right of this column, and those that are always red correspond to nodes to the left of this column. We note that the spatial mean cell density during an oscillation, plotted in Figure 6b, has many turning points despite being simply periodic. For this particular parameter set, changing the horizontal boundary conditions from no-flux on the fluid and cells to periodic conditions eliminated the oscillations (and solutions tended to a steady state). Similarly, the 1-D model did not display oscillating behaviour for these parameters.

In Figure 7a we demonstrate vertically asymmetric oscillations in an $n=25$ lattice. We plot the cell density at every node and demonstrate a similar oscillation as in Figure 5. We also show the effects of changing the nature of the boundary conditions, and the dimension of the domain, on the behaviour shown in Figure 7a. We show plots of cell density for periodic horizontal boundary conditions for both the fluid and the cell density in Figure 7b, and for the 1-D model in Figure 7c. The periodic conditions in Figure 7b change how the oscillation interacts with the boundary, and hence the bands do not display the same structure as in Figure 7a. The 1-D model in Figure 7c does not oscillate and instead reaches a steady state. The periods of the vertically asymmetric oscillations in Figures 5 and 7 are much larger than the vertically symmetric oscillations in Figure 3. In addition to the absence of asymmetric oscillations in the 1-D models, this change in oscillation frequency is how we differentiate these two kinds of oscillatory behaviours.

We also observe combinations of both kinds of oscillations. In Figure 8 we compare 1-D and 2-D simulations of a lattice of size $n=25$ with $p_{l}=750$ for different values of $\delta$. In 8b we observe asymmetric oscillations in a 2-D lattice in combination with smaller amplitude and higher frequency oscillations for some nodes. In 8d we instead observe fast oscillations with a slow variation indicative of two kinds of oscillatory behaviours. Figures 8a and 8c are the 1-D versions of Figures 8b and 8d respectively, and indicate the strength of the underlying pulsing oscillation in the 2-D simulations. Figure 8c shows more nodes oscillating with larger amplitudes compared to Figure 8a, which suggests why the oscillations in Figure 8d are qualitatively different from those in Figure 8b. Note that the frequencies of the fast oscillations in the 2-D simulations correspond to the frequencies of the 1-D oscillations, and hence the underlying vertically symmetric behaviour. We also find parameters where a pulsing oscillation in 1-D exists along with a vertically asymmetric steady state in the 2-D model, leading to a vertically asymmetric oscillation that is at the frequency of the vertically symmetric 1-D oscillation, but it is not vertically symmetric due to the underlying cell density distribution (i.e. the mean distribution is not vertically symmetric). In Section 5 we explore the parameter space in the 1-D and 2-D models, and demonstrate that regions with large oscillation amplitudes at some nodes in the 1-D model will correspond to pulsing oscillations in the 2-D model, so that these more exotic behaviours exist only when the 1-D model has smaller amplitude oscillations.

Some simulations exhibit excitable dynamics. In Figure 9a we first show the evolution of cell density for a small ($n=10$) lattice to steady state. Then in Figures 9b through 9d we add increasing perturbations to this steady state cell density at a single node at the top-right of the lattice. We observe small excursions from the equilibrium for very small perturbations (Figure 9b), larger excursions from equilibrium for larger perturbations (Figure 9c), and finally a sufficiently large perturbation induces a pulsing oscillation in Figure 9d. This excitability is due to the nonlocal nature of the quasi-static fluid coupling, in that an increase in cell density at one node in the lattice affects the overall flow through the scaffold instantaneously and hence the pressure at every other node. This is particularly true for nodes that have high cell density as the nonlinearity in Equations (1) and (3) amplifies the effect of large cell densities. Note in particular that we only perturb one node in Figures 9b-9d, yet several nodes grow or die rapidly at $t=20$ due to this nonlocal effect. In Section 5, we use numerical continuation to show the coexistence of locally stable steady states and pulsing oscillations for large regions of parameter space. A sufficiently large perturbation can move the state of the system between basins of attraction of these different long-time behaviours.

The smoothing parameter $g$ plays only a marginal role in these oscillations. The behaviours found in all of the simulations we performed are qualitatively and quantitatively the same in the limit of $g\to\infty$ (we replaced the functions $F_{1}$ and $F_{2}$ in Equations (6)-(7) with Heaviside step-functions, and obtained quantitatively comparable results). In the other direction, we set $g=1$ and repeated the simulation in Figure 3a, and the results displayed only small differences in the initial transient leading up to the oscillation, and otherwise the same qualitative behaviour. In contrast, the oscillations in Figures 5 and 7 were not present if $g=1$ for any parameters simulated, and instead steady state behaviour was obtained. This indicates that these vertically asymmetric oscillations are related to sharp transitions between cell growth and death due to pressure.

Aside from differences in the constitutive relations, the behaviour of the PDE model and larger lattices is similar, whereas the asymmetric solutions in smaller lattices do not have analogues in the PDE model. This is expected due to the similarities between the lattice equations and a method-of-lines discretization of the PDE. In this Section we have illustrated typical solution behaviours displayed by lattice and PDE models. In Section 5.2 we use simulations from a large subset of the feasible parameter space ($n$, $\delta$, $p_{l}$, and $p_{c}$) to understand where in parameter space these behaviours arise. In particular, we provide evidence that smaller lattices exhibit more complicated dynamics. Before demonstrating these results, we describe some asymptotic approaches that allow us to restrict our numerical efforts to specific parameter regimes.

4 Asymptotic solutions in the large diffusion limit

For larger values of the diffusion parameter $\delta$, and some values of $p_{l}$, the numerical solutions approach an equilibrium with little variation between nodal values (see Figures 2a and 2e where $\delta=10^{-1}$ is sufficient to observe very little spatial variation). We pursue an asymptotic solution for both models in the limit where $1/\delta\ll 1$. This asymptotic solution allows us to quantitatively predict the solution structure for large values of $\delta$. In particular, large-$\delta$ values give rise to simple steady state solutions, so that when we exhaustively explore a subset of the parameter space in Section 5.2, we only consider small values of $\delta$, namely $\delta\leq 10^{-1}$, in order to capture oscillations.

To make analytical progress we simplify the pressure forcing terms throughout the rest of this Section. Letting $g\to\infty$ our pressure functions become $F_{1}(p_{i})=H(p_{l}-p_{i})$ and $F_{2}(p_{i})=H(p_{i}-p_{l})$, where $H$ is the Heaviside step function. For convenience, we take $H$ to be continuous on the right, so that $\lim_{{\varepsilon}\to 0^{+}}H({\varepsilon})=H(0)\equiv 1$.

4.1 Lattice asymptotics

We first consider solutions in the large diffusion limit for the lattice. We rewrite Equation (2) as,

[TABLE]

We expand our variables as $N_{i}=N_{0,i}+\delta^{-1}N_{1,i}+O(\delta^{-2})$ and $p_{i}=p_{0,i}+\delta^{-1}p_{1,i}+O(\delta^{-2})$. Substituting these expansions into Equation (35) and formally equating powers of $\delta$, the leading-order equation for the cell density, $N_{0,i}$, is

[TABLE]

The graph Laplacian $L$ has one zero eigenvalue with the corresponding eigenvector $v=(1,1,,\dots,1)^{T}$ Newman [2010]. Therefore the leading-order cell densities are all equal; that is $N_{0,i}=N_{0}$ for some $N_{0}$, for all $i$.

To compute $N_{0}$ we consider the next-order problem,

[TABLE]

Equations (1) do not explicitly depend on $\delta$ and so the pressures $p_{0,i}$ satisfy Equations (1) with the uniform cell density $N_{0}$. As a spatially-constant cell density $N_{0}$ is vertically symmetric, we can use Equation (23) to determine the pressures $p_{0,i}$. In order to eliminate the $N_{1,j}$ terms, we sum Equation (37) over $i$ and note that $\sum_{i=1}^{n^{2}}\sum_{i=1}^{n^{2}}L_{ij}N_{1,j}=0$, to find

[TABLE]

where we have relabelled the pressures to reflect the 1-D symmetry, and divided both sides by $n$.

The pressures $\hat{p}_{0,i}$ are monotonically decreasing in $i$ with $\hat{p}_{0,n}=1/n$. Excluding cases corresponding to $p_{l}$ outside of the ranges given by (33), there exists an $m=m(N_{0})$ such that for all $j\leq m$ and $k>m$, $\hat{p}_{0,j}\geq p_{l}$ and $\hat{p}_{0,k}<p_{l}$. Thinking of $m$ as an integer, we have that solutions of Equation (38) satisfy,

[TABLE]

We now approximate $m(N_{0})$ as a real variable, rather than an integer, to make analytical progress and avoid many tedious details. We will justify this approximation later by observing how well the solution matches numerical solutions of the full problem. By setting the left hand side of Equation (23) equal to $p_{l}$, we find $m(N_{0})$ to be

[TABLE]

where $[f]^{+}=\max(f,0)$ is the positive part of the function $f$. It is necessary to define $m$ using the notation $[f]^{+}$ as for small $N_{0}$, corresponding to $p_{i}<p_{l}$ for all $i$, the value of $m(N_{0})$ is negative. Equations (39) and (40) represent a single scalar first-order ordinary differential equation, compared with the original $n^{2}$ differential equations (2) and $n^{2}$ algebraic equations (1), although our scalar equation is not smooth (e.g. not differentiable) due to the function $m(N_{0})$.

We now look for steady state solutions, $N_{0}^{*}\geq 0$, to Equation (39). If $m(N_{0}^{*})=0$, then $N_{0}^{*}=0$ or $N_{0}^{*}=1$ are the only solutions to Equation (39), and we will show that they cannot be stable steady states within the parameters we consider. Assuming that $m(N_{0}^{*})>0$, we substitute (40) into Equation (39), and find that $N_{0}^{*}$ satisfies the quintic polynomial

[TABLE]

We will first argue that any root of Equation (41) must be unique within the interval $(0,1)$, and then show the existence and global stability of $N_{0}^{*}\in(0,1)$.

Let $\tilde{N}=R(N_{0}^{*},N_{0}^{*})=1-\rho N_{0}^{*}$. Since $0<\rho<1$, all roots $N_{0}^{*}\in(0,1)$ of Equation (41) must be associated to roots $\tilde{N}\in(0,1-\rho)$ of

[TABLE]

and vice versa. Equation (42) has only one sign change in its coefficients (as $p_{l}>1/n$ by (33)). By Descartes’ rule of signs, Equation (42) has at most one positive root, implying that there cannot be more than one root in the interval $0\leq\tilde{N}\leq 1-\rho$. Hence any root of equation (41) in $(0,1)$ must be unique.

We now argue by monotonicity that this root exists and is globally attractive for initial data in the full interval $[0,1]$. The function $m(N_{0})$ given by (40) is monotonically increasing in $N_{0}$, and for $1+1/n<p_{l}<(1+(n-1)(1-\rho)^{-4})/n$ we have that $m(0)=0$ and $m(1)>0$. So there is a unique value $\hat{N_{0}}$ such that $m(\hat{N_{0}})=0$ and for all $N>\hat{N_{0}}$, $m(N)>0$. From Equation (40) we compute that

[TABLE]

which is within the interval $[0,1]$.

By monotonicity, we have that $m(N_{0})\geq 0$ for all $N_{0}\in[\hat{N_{0}},1]$. Since $m(0)=0$ and $m(1)>1$, we have that $N_{0}=0$ is an unstable steady state and $N_{0}=1$ is not a steady state solution to Equation (39) for these parameters. From Equation (39) we have that $N_{0}^{\prime}(t)<0$ when $N_{0}=1$ and $N_{0}^{\prime}(t)>0$ when $N_{0}=0$. This implies that continuous solutions to Equation (39) do not leave the interval $[0,1]$, which is consistent with the boundedness of solutions to the full system shown in Section 2.3. We note that for $m(N_{0})=0$, the right side of (39) is positive, so $N_{0}(t)$ grows in time and approaches $\hat{N_{0}}$. Similarly, for $N_{0}=\hat{N_{0}}$, the right side of (39) is positive as $1-\hat{N_{0}}>0$. So the region $(\hat{N_{0}},1)$ absorbs solutions with initial data $N_{0}(0)\in(0,1]$, and hence we must have that there exists a steady state solution $N_{0}^{*}\in(\hat{N_{0}},1)$ to (39). Both existence and global attractivity of this steady state follows from the fact that Equation (39) is just a scalar first-order ODE within the region $N_{0}^{*}\in(\hat{N_{0}},1)$.

We can solve equation (39) by considering the cases when $N_{0}(0)<\hat{N_{0}}$ and $N_{0}(0)>\hat{N_{0}}$ separately. We first assume that $N_{0}(0)\in(0,\hat{N_{0}})$. By the monotonicity of $m$, and the fact that $N^{\prime}(t)>0$ as long as $m(N)=0$, there exists some $t_{s}$ such that $N_{0}(t_{s})=\hat{N_{0}}$. So substituting (40) into (39) we have

[TABLE]

The solution to equation (44) for $t\leq t_{s}$ is a logistic function,

[TABLE]

which we use to compute $t_{s}$. We set $N_{0}(t_{s})=\hat{N_{0}}$ and use equations (43) and (45) to find

[TABLE]

We consider equation (44) for $t>t_{s}$. This equation has a sixth-order polynomial nonlinearity, and we are not aware of a method to explicitly integrate it. We know by the boundedness of $N_{0}(t)$, however, that $N_{0}(t_{s})$ is nearby in a nondimensional sense to our equilibrium value $N_{0}^{*}$ for any initial mean cell density $N_{0}(0)$. So we expect approximately exponential convergence to $N_{0}^{*}$ for $t>t_{s}$. We substitute the ansatz

[TABLE]

for $\left|A\right|\ll 1$ into equation (44), and take the Taylor series of the right hand side about $A=0$ to find

[TABLE]

The $O(1)$ term on the right hand side of equation (48) is zero, as $N_{0}^{*}$ is the solution to equation (41), which is obtained numerically. Comparing $O(A)$ terms we have

[TABLE]

In order to make the solution (47) continuous with our solution for $t\leq t_{s}$, we take $A=N_{0}(t_{s})-N_{0}^{*}$ (which for our dynamics satisfies $|N_{0}(t_{s})-N_{0}^{*}|\leq 1$) where $N_{0}(t_{s})$ is the solution from equation (45) evaluated at $t=t_{s}$. Our asymptotic solution can then be written as

[TABLE]

where $t_{s}$ is defined by equation (46) and $r$ by equation (49).

Next we consider the case that $N_{0}(0)>\hat{N_{0}}$. We have that $m(N_{0}(0))>0$, and so we can use the method above to find an exponential solution of the form given by equation (47) as before. The value of $A$ is again determined by forcing the solution to continuously approach the initial data, which gives that $A=N_{0}(0)-N_{0}^{*}$. So the solution in this case is

[TABLE]

with $r$ given by equation (49).

The solution given by equation (50) describes a period of logistic cell growth, followed by an exponentially fast convergence to a steady state, and depends on the initial mean cell density $N_{0}(0)$, $p_{l},\rho$ and $n$. If the initial mean cell density is high enough to immediately cause cells to die due to pressure, there is no logistic growth phase and Equation (51) describes the exponential convergence to a steady state. These analytical solutions also match numerical solutions of (39) extremely well; see Figure 10. In the worst case of $N_{0}(0)=0.9$ and $p_{l}=10$, where there is no logistic growth, we note that the analytical solution is still quite close to the numerical solution despite the value of $A=0.45$, which is not asymptotically small. We suspect that this close agreement is due in part to solutions being bounded in $[0,1]$, so that we always have $\left|A\right|<1$, along with the existence of a unique stable steady state to equation (39).

We compare the asymptotic solution (50) for the cell density, $N_{0}$, against the (spatial) mean cell density for the full two-dimensional problem (1)-(7) for various values of $\delta$ in Figure 11. For small values of $\delta$, the oscillations in nodal cell densities have an overall effect on the spatial mean cell density such that our asymptotic solution over-predicts the mean cell density of oscillatory solutions. Increasing $\delta$, we see a relatively fast convergence to the asymptotic solution for the spatial mean cell density $N_{0}$. Comparable results hold for different values of $p_{l}$ and $n$, and for the solution given by equation (51).

In Figure 12 we plot values of (spatial) mean cell densities for simulated solutions of equations (1)-(7) for a lattice with $n=10$, along with corresponding $N_{0}^{*}$ computed from Equation (41). In both cases we plot the steady state behaviour, so if the simulations of the full problem were oscillatory, we averaged over $t=10$ units of time. The large $\delta$ steady state $N_{0}^{*}$ accurately predicts the long-time spatial mean cell density for $\delta\approx 1$ for all values of $p_{l}$, and for smaller values of $\delta$ for larger values of $p_{l}$. For smaller values of $\delta$, the variation from the asymptotic solution is partly due to significant spatial variation in the cell density $N_{i}$. The jagged regions in the plots are due to changes in the number of grid points that are growing or dying, rather than coarse plotting. The convergence in $\delta$ being faster for larger values of $p_{l}$ explains why in Figure 2a the cell density is very close to uniform, whereas in Figure 2b there is a substantial deviation from uniformity.

4.2 Continuum asymptotics

Analytical solutions for the PDE model can be obtained in exactly the same way as in Section 4.1, so we just mention key steps here. We expand our variables as $N=N_{0}+\delta^{-1}N_{1}+O(\delta^{-2})$ and $p=p_{0}+\delta^{-1}p_{1}+O(\delta^{-2})$. Substituting these into Equation (11), we find the leading-order equation for cell density,

[TABLE]

which, along with the boundary condition in equations (19), shows that $N_{0}$ is constant throughout the domain (e.g. by the Maximum Principle). To find $N_{0}$ we consider the next order problem,

[TABLE]

We now eliminate the dependence on $N_{1}$ using the divergence theorem. Integrating equation (53) over the domain, along with the boundary conditions in (19), we find

[TABLE]

where $p_{0}(x)$ is a function of $x$ only due to the leading-order cell density $N_{0}$ being constant.

From equation (29) we have that $p_{0}(x)=(1-x)/k(N_{0})$ which is monotonically decreasing in $x$. Given the bound from (34), which implies that we will see neither uniform logistic growth or uniform cell death, there must exist a point $x^{*}(N_{0})$ such that for all $x\leq x^{*}(N_{0})$, $p_{0}(x)\geq p_{c}$, and for all $x>x^{*}(N_{0})$, $p_{0}(x)<p_{c}$. From Equation (54) we then have,

[TABLE]

where $x^{*}(N_{0})=[1-p_{c}k(N_{0})]^{+}=[1-p_{c}(1-\rho N_{0})^{3}]^{+}$, where again the positive part must be taken to avoid a negative value of $x^{*}(N_{0})$ for small $N_{0}$. We have reduced a system of two partial differential equations in equations (9)-(11) to equation (55) which is an ordinary differential equation, although $x^{*}(N_{0})$ is not everywhere-differentiable. At steady state, the solution $N_{0}^{*}$ must satisfy,

[TABLE]

where we are neglecting the unstable steady state $N_{0}^{*}=0$. This is a quartic equation for $N_{0}^{*}$ which can be solved numerically. Up to constitutive differences giving different exponents between $k(N_{0})$ and $R^{-4}(N_{0})$, it can also be seen as the limit of equation (41) as $n\to\infty$. As in equation (41), equation (56) can be shown to have at most one solution, $\tilde{N}$, in the interval $(0,1)$ which is globally attractive.

We can follow the same procedure that was employed in Section 4.1 to derive Equations (50) and (51) to solve equation (55). We again consider a period of logistic growth (assuming the initial cell density is not too large), wherein the solution can be found analytically. As the cell density becomes large enough to induce non-uniform death (e.g. $x^{*}(N_{0}>0)$), we exploit the fact that the difference between $N_{0}(t)$ and the equilibrium solution given by Equation (56) is small to find the solution for all time $t$. We define $\hat{N_{0}}$ to be the value of cell density such that $x^{*}(\hat{N_{0}})=0$ and for all $N_{0}>\hat{N_{0}}$, $x^{*}(N_{0})>0$. If we assume that $N_{0}(0)<\hat{N_{0}}$, then the solution takes the form

[TABLE]

where $N_{0}(0)$ is the initial mean cell density, $N_{0}^{*}$ is the unique root of equation (56), the switching time $t_{c}$ is given by

[TABLE]

and $r_{c}$ is given by

[TABLE]

Similarly, if we assume that the initial cell density satisfies $N_{0}(0)>\hat{N_{0}}$, then the solution is

[TABLE]

Solutions (57) and (60) are analogous to the lattice solutions (50) and (51). In Figure 13 we compare the analytical and numerical solutions to equation (55). The worst case behaviour, shown in red, corresponds to the asymptotic parameter $N_{0}^{*}-N_{0}(t_{c})=0.5$, and it is not surprising that there is a transient discrepancy between the numerical and analytical solutions in this case.

We compare this solution for the cell density $N_{0}$ against the spatial mean cell density for the full 2-D problem (9)-(19) for various values of $\delta$ in Figure 14. For small values of $\delta$, as in the lattice case, there is some discrepancy between the numerical and asymptotic solutions due to spatial structure. As we increase $\delta$, however, we see a relatively fast convergence to the asymptotic solution for $N_{0}$. Comparable results hold for different values of $p_{c}$. We note that the asymptotic solution converges more quickly to the equilibrium value $N_{0}^{*}$ than the numerical solution of Equation (55).

In this section, we have shown that there exist parameter regimes with relatively simple dynamics for large values of $\delta$, or moderate $\delta$ and large values of $p_{l}$ or $p_{c}$ (see Figure 12; the PDE behaviour is qualitatively similar). In these regimes, the transient and long-time dynamics of both the lattice and PDE models can be captured by simple scalar ODEs each admitting a unique attracting steady state (Equations (44) and (55) respectively). We anticipate that the oscillatory solutions observed in Section 3 can be confined to regions of the parameter space where either $\delta$ is small, or the threshold $p_{l}$ or $p_{c}$ is small. We also note that other simulations show that the dimensionality (e.g. 1-D or 2-D) or the nature of the boundary conditions of the model (e.g. periodic) does not play a role in solution behaviour within this asymptotic regime (large $\delta$), whereas these model details do play a role in the nature of oscillations and steady states observed outside of this asymptotic regime (see Figure 7).

5 Bifurcations in small and large lattices

We now consider the behaviours observed in Section 3 from the perspective of how solution behaviours change as parameters are varied. We begin by demonstrating the existence of Hopf and pitchfork bifurcations from steady states in an $n=4$ lattice as $\delta$ varies, and describe how the size of the lattice affects the existence of these bifurcations. Finally we classify behaviours observed in simulations throughout a bounded subset of the parameter space of $\delta$ and $p_{l}$ (or $p_{c}$) for larger lattices and the PDE.

5.1 Existence of symmetry-breaking bifurcations

We examine a low-dimensional ($n=4$) lattice model, and numerically continue steady state solutions over a range of $\delta$ with $p_{l}=15$. We track the local stability of these steady state solutions to small amplitude perturbations by evaluating the Jacobian of the system and its eigenvalues. We use natural parameter continuation and manually switch between solution branches by continuing solutions in $\delta$ from any stable equilibria found by solving Equation (2). That is, we take our equilibrium solution for a specific value of $\delta$, $N_{i}^{*}(\delta)$, then solve the algebraic steady state Equation (2) using the Matlab function ‘fsolve’ with $N_{i}^{*}(\delta)$ as an initial guess for the solution at a new value of the parameter $\delta^{\prime}=\delta\pm\epsilon$ for some small $\epsilon$. We choose $\epsilon=10^{-4}$, and continue solutions forward and backward for all $\delta\in[0,2]$.

Equation (1) can be viewed as a weighted graph Laplace equation, and can be inverted numerically to find the pressures as functions of the cell density. We consider the $n^{2}\times n^{2}$ matrix corresponding to the Jacobian of equation (2) with $p_{i}$ a function of $N_{1},\dots,N_{n^{2}}$ for all $i$. We denote the $ith$ eigenvalue of the Jacobian evaluated along a branch of steady states by $\sigma_{i}$. See Seydel [2009] for discussion of continuation procedures in general.

In Figure 15, we plot the value of the cell density at the bottom-left node, $N_{1}$, corresponding to different steady state solution branches, for varying $\delta$ found from the above procedure. The solid lines denote steady states with all eigenvalues having negative real part, $\max_{i}\{\Re(\sigma_{i})\}<0$, and hence are locally linearly stable. The dashed lines have at least one eigenvalue with positive real part, and so are linearly unstable. The black line corresponds to a vertically symmetric equilibrium (which remains vertically symmetric for all $\delta$), and the red and blue lines to vertically asymmetric equilibria (see the insets of Figure 15 which show cell density distributions for these equilibria at $\delta=0.5$).

Around $\delta\approx 0.08$ (point A), the vertically symmetric steady state loses stability. Beyond point A there is an oscillating solution that attracts solutions perturbed from the now unstable vertically symmetric solution. This oscillating solution is of the vertically asymmetric kind plotted in Figure 5. We plot the envelope of the oscillation that $N_{1}$ undergoes during a period. We plot the two eigenvalues with largest real part in Figure 16, and show that at the point where the vertically symmetric steady state loses stability, the eigenvalues cross the imaginary axis in a conjugate pair (e.g. $\Im(\sigma_{i})\neq 0$ for the two complex conjugate $\sigma_{i}$ with $\Re(\sigma_{i})=0$). Since the oscillatory solution immediately beyond this loss of stability is numerically stable and its amplitude is small for $\delta$ near the bifurcation point, this indicates a supercritical Hopf bifurcation giving rise to the stable limit cycle that we observe for $\delta>0.08$ in Figure 15.

For $\delta\approx 1.5$ (point D), two solution branches emerge from the vertically symmetric branch as it loses stability (for decreasing $\delta$). These vertically asymmetric solutions are reflected copies of one another, as is typical in a pitchfork bifurcation. We plot the eigenvalues of these continued branches in Figure 17 showing evidence that this is a pitchfork bifurcation. The eigenvalues touch the imaginary axis with zero imaginary part at D in Figure 17 then return to the left hand side of the plane, whereas in Figure 16 the eigenvalues cross the imaginary axis as the vertically symmetric solution loses stability. The asymmetric branches both lose stability at $\delta\approx 0.46$ (point C), again satisfying the criteria for a Hopf bifurcation, as the eigenvalues with the largest real part again cross the imaginary axis in Figure 17 (but here for decreasing $\delta$). In this case, however, we do not find a small amplitude stable limit cycle, and instead solutions are attracted to the limit cycle previously generated, which has a large amplitude at this point. This shows that there is a subcritical Hopf bifurcation for both of the asymmetric equilibria at C, which leads to the creation of an unstable limit cycle.

For larger values of diffusion the limit cycle is no longer observable and we suspect it becomes unstable around $\delta\approx 0.47$ (point C’). We know that an unstable limit cycle exists to the right of C, between C and C’, and that this may lead to the loss of stability for the stable limit cycle at C’, but we do not pursue these claims here. Lastly, we see that around $\delta\approx 0.4$ (point B) the three steady state solution branches merge, and the eigenvalues shown in Figure 17 have the same structure for the rest of the continued solution. Note that near B all branches have large real eigenvalues ($\sigma_{i}\approx 20$), and so there is no exchange of stability. While we have only shown plots for $N_{1}$, plots for all $N_{i}$ show that these behaviours are qualitatively consistent across all states, and hence these branches, their stability, and the merging of branches are faithfully captured by Figure 15.

The results in Figures 15-17 suggest the existence of Hopf and pitchfork bifurcations, alongside multiple steady states, and oscillatory behaviour in the 2-D lattice with $n=4$. We did not observe any pulsing oscillations of the type shown in Figure 4. We briefly summarize results for other values of $n$. For $n=2$ and $n=3$ we only observed steady states for a large region of parameter space. Similar behaviours to $n=4$ (vertically asymmetric and symmetric steady states and asymmetric oscillations) were observed for $n=5,6$. For $n=7$, the 2-D lattice had both vertically symmetric pulsing oscillations, as in Figure 4 and vertically asymmetric oscillations, in addition to steady state behaviour. For the 1-D model given by Equations (24), we only observed steady state behaviour for $n\leq 6$, and oscillations and steady states for $n\geq 7$.

For parameter values where a vertically symmetric oscillation exists but a vertically asymmetric oscillation does not, the vertically symmetric steady state remains locally stable. The 1-D model exhibits a locally stable steady state solution for every parameter combination we simulated. Together, this suggests that it is only these asymmetric oscillations or steady states that can induce instability in vertically symmetric 2-D steady states, as these vertically symmetric steady states corresponds to steady states in the 1-D model which are always (locally) stable. The multiple stability of steady states and oscillations (e.g. between point C and C’ in Figure 15) observed also provides an explanation for the excitable behaviour in Figure 9, as the perturbations need to be sufficiently large to move out of the basin of attraction of the locally stable steady state.

These results suggest that the vertically asymmetric oscillations in the smaller lattices are due to local symmetry-breaking bifurcations of Hopf type. For general differential-algebraic systems, codimension-1 bifurcations can also be due to a singularity in the algebraic subsystem, or due to saddle-node bifurcations Venkatasubramanian et al. [1995]. By the regularity of the graph Laplacian in (1), we conjecture that no singularity-induced bifurcations can occur. While we are unable to rule out saddle-node bifurcations, we do not observe them in any of our continuation studies. The vertically symmetric pulsing oscillations are not, as far as we can detect numerically, created via a Hopf or other local bifurcation from a steady state.

5.2 Classification of the parameter space

Even for small lattices, a complete classification of the parameter space is not tractable. Instead we now consider lattices of varying size and classify behaviours broadly as oscillatory or steady-state, and record properties of solutions. Namely, we record the amplitude and frequency of oscillations. The results from Section 4 suggest that for large values of diffusion the cell density is approximated accurately by a spatially constant solution. For this reason we choose to do parameter sweeps between $\delta=0.001$ and $\delta=0.1$, and variations in the threshold parameters within the bounds given by (33)-(34). We compute bifurcation plots in $(\delta,p_{l})$ and $(\delta,p_{c})$ in the following way. We first discretize the parameter space, and numerically simulate the governing equations at each discrete point for a sufficiently long period of time (we set this time to be $t=40$, and use identical initial conditions described at the start of Section 3). We then truncate our simulation to analyze only the last $10$ units of time. For each node, or in the case of the PDE for each interpolated discretized element, we compute the largest and smallest values this node takes in the truncated time series to compute a nodal oscillation amplitude.

In Figure 18 we plot the maximal (over all nodes) amplitude of these nodal oscillations in the lattice, and in Figure 19 we plot the frequency of nodal oscillations computed using the Fast-Fourier Transform of these truncated time series. We plot corresponding PDE bifurcation diagrams in Figure 20. The larger lattices in Figure 18 and the PDE in Figure 20 appear to have an oscillatory region for small $\delta$ or small $p_{l}$. The 2-D lattices have more complicated regions not present in the 1-D models where vertically asymmetric oscillations are present (compare the left and right columns of Figure 18, especially for $n\geq 50$). The magnitude of the nodal oscillations away from the pulsing regime (bottom-right region of each plot in the right column of Figure 18) appears to decrease as a function of $n$, whereas the pulsing oscillations maintain a comparable amplitude for $n\geq 25$ (not shown). Regions with large (nodal) amplitudes in the 1-D model correspond to pulsing oscillations in the 2-D model. The mixed oscillation (vertically symmetric and asymmetric) behaviours shown in Figure 8 exist only when the 1-D model has smaller amplitude oscillations, and hence along the boundaries of the oscillatory regions in the plots along the left of Figure 18. Note too that for $n=10$ both the 1-D and 2-D oscillation regions are much more disconnected than the corresponding regions for larger lattices (this disconnected bifurcation structure also existed for smaller lattices). For $n\geq 25$ the region where oscillatory behaviour is observed appears to become simply connected in the 1-D models, but it is unclear that this happens in the 2-D case for large $n$.

The frequency of oscillations away from the small $\delta$ or $p_{l}$ regime increases with $\delta$ (see Figure 19). The oscillation frequency is substantially larger for oscillations present in the 2-D lattices that are not present in their 1-D counterparts; compare the colour ranges in Figure 19. The PDE behaviours shown in Figure 20 have no appreciable difference between the 1-D and 2-D model behaviours. This suggests that the symmetry-breaking oscillations of the smaller lattices may disappear in the limit of $n\to\infty$.

We also consider the amplitudes of oscillations in the spatial mean cell density over the entire lattice in Figure 21, in comparison to the maximal nodal oscillations plotted in Figure 18. We see no qualitative difference between the solutions of the 1-D and 2-D models, suggesting that for large $n$, the asymmetric oscillations become less influential on the spatial mean cell density. This is also consistent with the maximal amplitudes of the nodal oscillations decaying with $n$.

6 Discussion

In this paper we have presented and analyzed a pair of new models based on those proposed in Krause et al. [2017] that describe an active porous medium with interactions between cell growth and fluid flow. We have demonstrated several notable differences between the lattice models of various sizes, as well as between the lattice models and the spatially continuous model.

In Section 2 we presented the lattice and continuum models, and demonstrated a reduction to a 1-D model available for vertically symmetric cell density distributions. In Section 3, we solved the (1-D and 2-D) systems numerically to demonstrate the existence of steady state solutions and oscillations, either preserving or breaking vertical symmetry. Vertically symmetric solutions were robust to changing the boundary conditions and dimension of the models, as well as the size of the lattice above a minimum threshold. Vertically asymmetric solutions appeared in smaller lattices, and these were more sensitive to the nature of the boundary conditions (and did not exist in the 1-D models). In Section 4 we derived large $\delta$ asymptotic solutions for cell density in both the continuum and lattice models. In Section 5 we exhibited Hopf bifurcations from steady states to vertically asymmetric oscillations. The vertically symmetric oscillations, on the other hand, were not created from local bifurcations from steady states that we found numerically. The asymptotic solutions from Section 4 allowed us to confine parameter sweeps in Section 5 to ranges of the model parameters where variation in spatial and temporal structure may be found, and we characterized model behaviours across these parameter ranges. We found that there is a region of parameter space for small $\delta$ or small pressure thresholds where oscillations are observed. This region appears simply connected for large lattices, and for the continuum model, but for smaller lattices the parameter space consists of disconnected regions corresponding to oscillatory behaviour. Smaller lattices also exhibit both vertically symmetric and asymmetric oscillations.

Results in Sections 3-5 demonstrate an important role played by vertically symmetric steady states and oscillatory solutions to these models. Our numerical explorations show that in various limits, notably the large lattice limit of $n\to\infty$, these two solutions are the main stable long-time behaviours of the models, at least within the parameter regimes we explored. While asymmetric oscillations do exist in large lattices, they become almost indistinguishable from the vertically symmetric oscillations. For large lattices, quantities such as mean cell density oscillation amplitudes, shown in Figure 21, are completely indistinguishable. For smaller lattices we have shown the existence of Hopf and pitchfork bifurcations that break the symmetry of the vertically symmetric steady state solutions and give rise to asymmetric steady states, and asymmetric oscillating solutions. The phase space and the parameter space for moderate lattice sizes (e.g. $n\sim O(10)$) is difficult to analyze as it admits all of the various solutions, alongside what are complicated combinations of them (as demonstrated in Figure 8).

We can interpret our results in terms of solution symmetry and nonlocal reaction-diffusion mechanisms. Due to results from Casten & Holland [1978]; Matano et al. [1978], and others (see Smith [2008] for a more thorough review), it is known that local scalar reaction-diffusion equations with Neumann data on a bounded domain have relatively simple asymptotic behaviour - namely, all bounded solutions converge to spatially uniform equilibria. Equation (30) is a simplification of our model that still admits spatially-structured and oscillatory solutions despite being a scalar reaction-diffusion equation, and hence this nonlocality is a necessary component for both the spatially structured steady states and oscillatory behaviours. Simpler models of this kind of nonlocal reaction-diffusion have been considered in Billingham [2004]; Gourley et al. [2001]; Hamel & Ryzhik [2014]; Levchenko et al. [2016], and the references therein, where oscillations and non-uniform equilibria are found and analyzed. We note that these models are qualitatively similar to Equation (30) with different nonlinearities, both in the reaction kinetics and the integral kernel.

We conjecture that it is this inherently nonlocal structure, driven by the long-range quasi-static pressure forcing, that gives rise to the oscillatory behaviours we observe. For smaller lattices, the vertically asymmetric steady states and oscillations are related to the discrete structure of these finite-size lattices, as these are not present in the continuum analogue. We conjecture that this discrete structure leads to symmetry-breaking effects due, for example, to the anisotropic nature of discrete diffusion on a lattice Cardy [2012]. We conclude that the vertically asymmetric oscillations are then an interplay of this symmetry breaking with the nonlocal mechanism described before.

Our observations about vertically asymmetric oscillations are consistent with the oscillations observed in Krause et al. [2017], where a similar model captured cell death due to high values of shear stress. While the solution behaviours in Krause et al. [2017] were generally more complicated (e.g. lacked the spatial symmetries present here), the models there satisfied the same fluid equations and hence also had a nonlocal coupling between cell growth and fluid pressure as demonstrated explicitly in Sections 2.1 and 2.2. Oscillations in cell density were also observed in the lattice model of Krause et al. [2017], with the amplitude of oscillations decreasing with increasing $n$, but no oscillations were not observed in the spatially continuous analogue. By comparison, the spatially continuous pressure model in Section 2 did give rise to oscillations. The difference between these models may be due to a localization of the nonlocal pressure in the fluid shear stress, which is a function of gradients in pressure. This is in contrast to the models presented in Section 2, as the pressure itself is forcing the cell growth problem, and hence oscillations are observed in both lattice and continuum models. Overall, these results suggest that oscillations observed in Krause et al. [2017] are analogous to the vertically asymmetric oscillations in the pressure-forced model, as these are also a finite-size effect observed in smaller lattices.

In the context of bioactive porous media, these models demonstrate several interesting effects due to the finite pore network captured in our lattice model. In the limit of very large lattices, we suspect that spatially continuous models are good approximations to these kinds of pore networks, at least for the simple topology considered here Verwer & Sanz-Serna [1984]. For smaller networks, we have demonstrated nontrivial effects due to the discrete symmetries inherent in the finite-scale lattice that cannot, as far as we are aware, be captured in continuous analogues Cardy [2012]. We suspect that more realistic (and hence more complicated) pore networks may exhibit other behaviours due to the presence or lack of symmetries in the topology of the pore network. The tissue engineering literature, for instance, emphasizes macroscopic and microscopic properties of porous scaffolds, such as permeability and pore size respectively, but does not seem to consider as much the importance of the discrete structure of pore networks Hutmacher [2000]; Jones et al. [2008]; Wu et al. [2014]. Further explorations of network models of bioactive porous media, such as those studied here, will be both quantitatively useful for applications, and provide a number of interesting mathematical questions to pursue.

Acknowledgments

D.B. was partially funded by Engineering & Physical Sciences Research Council (EPSRC) Fellowship ref. EP/M002896/1. S.L.W. gratefully acknowledges funding in the form of a Royal Society Leverhulme Trust Senior Research Fellowship. In compliance with EPSRC’s open access initiative, the data in this paper is available from http://dx.doi.org/xxx/xxx.