Lattice and Continuum Modelling of a Bioactive Porous Tissue Scaffold

TL;DR

This paper compares lattice-based and continuum models of bioactive porous tissue scaffolds, highlighting the importance of explicitly modeling scaffold network structure for accurate predictions in tissue engineering.

Contribution

It introduces a lattice model that explicitly incorporates scaffold network structure and compares it with existing continuum models, revealing significant differences in predictions.

Findings

Lattice model captures spatial pattern formation and oscillations not seen in continuum models.

Explicit network modeling affects predictions of cell growth and fluid flow.

Differences depend on scaffold geometry and cell types.

Abstract

A contemporary procedure to grow artificial tissue is to seed cells onto a porous biomaterial scaffold and culture it within a perfusion bioreactor to facilitate the transport of nutrients to growing cells. Typical models of cell growth for tissue engineering applications make use of spatially homogeneous or spatially continuous equations to model cell growth, flow of culture medium, nutrient transport, and their interactions. The network structure of the physical porous scaffold is often incorporated through parameters in these models, either phenomenologically or through techniques like mathematical homogenization. We derive a model on a square grid lattice to demonstrate the importance of explicitly modelling the network structure of the porous scaffold, and compare results from this model with those from a modified continuum model from the literature. We capture two-way coupling between cell growth and fluid flow by allowing cells to block pores, and by allowing the shear stress of the fluid to affect cell growth and death. We explore a range of parameters for both models, and demonstrate quantitative and qualitative differences between predictions from each of these approaches, including spatial pattern formation and local oscillations in cell density present only in the lattice model. These differences suggest that for some parameter regimes, corresponding to specific cell types and scaffold geometries, the lattice model gives qualitatively different model predictions than typical continuum models. Our results inform model selection for bioactive porous tissue scaffolds, aiding in the development of successful tissue engineering experiments and eventually clinically successful technologies.

Tables

Table 1. Table 1 : Parameters used in the models. 1 Vunjak-Novakovic et al., ( 1998 ) , 2 Vafai, ( 2010 ) , 3 Shakeel et al., ( 2013 ) , 4 McCoy et al., ( 2012 ) , 5 Truscello et al., ( 2012 ) , 6 Glowacki et al., ( 1998 ) , 7 Whittaker et al., ( 2009 ) , 8 Koponen et al., ( 1996 ) , 9 Kou et al., ( 2012 ) , .

Parameter	Expression	Values
Cell volume1	$\nu$	$4.2\times {10}^3-1.7\times {10}^4\mu m^3/\text{cell}$
Maximum cell density 1 2 3	$N_c$	${10}^{13}-{10}^{14}\text{cell}/m^3$
Scaffold length1	$L$	$2-10mm$
Pore length scale4	$l,R_0$	$85-325\mu m$
Cell proliferation timescale2 3	$1/\beta$	$6.59\times {10}^{4}s$
Viscosity1	$\mu$	$0.7\mathrm{mPs}$
Cell-free scaffold porosity3 5	${\varphi }_0$	$40$-$97$%
Pump flow rate4 6	$Q_c$	$5.2-6.24{\mathrm{mm}}^2/s$
Cell diffusion rate2 3	$D_n$	${10}^{-12}-{10}^{-11}{m}^2/s$
Tortuosity7 8 9	$\tau$	1.15
Ratio of timescales	$\delta =D_n/L^2\beta$	${10}^{-4}-{10}^1$
Maximum cell volume fraction	$\rho =\nu N_c$	${10}^{-2}-1$
Pores per lattice side	$n=L/l$	$10-{10}^3$

Equations

$${\mathbf{u}}^{\ast }=-\frac{k^{\ast }(x^{\ast },y^{\ast },N^{\ast })}{\mu }{\nabla }^{\ast }{p}^{\ast }\text{and }{\nabla }^{\ast }\cdot {\mathbf{u}}^{\ast }=0,$$

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

$$p^{\ast }=p_0\text{at}{x}^{\ast }=0,0\le y^{\ast }\le L,\text{and}{p}^{\ast }=p_1\text{at}{x}^{\ast }=L,0\le y^{\ast }\le L,$$

$$Q_0^{\ast }={\int }_0^L-\frac{k^{\ast }(0,y^{\ast },N^{\ast })}{\mu }\frac{\partial p^{\ast }}{\partial x^{\ast }}d{y}^{\ast }.$$

$${{\mathbf{u}}_{\mathbf{r}}}^{\ast }=\frac{Q_c}{Q_0^{\ast }}{\mathbf{u}}^{\ast }.$$

$$\frac{\partial N^{\ast }}{\partial t^{\ast }}=\beta \left(F_1({\sigma }^{\ast })N^{\ast }\left(1-\frac{N^{\ast }}{N_c}\right)-F_2({\sigma }^{\ast })N^{\ast }\right)+{\nabla }^{\ast }\cdot (D^{\ast }(N^{\ast }){\nabla }^{\ast }{N}^{\ast }),$$

$$\mathbf{n}\cdot \nabla N^{\ast }=0\text{for }\mathbf{x}\in \partial [0,L{]}^2,$$

$$N^{\ast }(x^{\ast },y^{\ast },0)=N_0^{\ast }(x^{\ast },y^{\ast }),$$

$$\hat{\varphi }(N^{\ast })={\varphi }_0(1-\nu N^{\ast }),$$

$$k^{\ast }(N^{\ast })=k_0\hat{\varphi }(N^{\ast }{)}^3,$$

$${\sigma }^{\ast }=\frac{4\mu \tau }{R_0}\frac{\Vert {{\mathbf{u}}_{\mathbf{r}}}^{\ast }\Vert }{\hat{\varphi }(N^{\ast })},$$

$$F_1({\sigma }^{\ast })=1-\left(\frac{1}{2}\right)(\tanh [g({\sigma }^{\ast }-{\sigma }_t)]+1),\text{and }{F}_2({\sigma }^{\ast })=\left(\frac{1}{2}\right)(\tanh [g({\sigma }^{\ast }-{\sigma }_t)]+1),$$

$$(x^{\ast },y^{\ast })=L(x,y),{\nabla }^{\ast }=\frac{1}{L}\nabla ,(N^{\ast },N_0^{\ast })=N_c(N,N_0),\hat{\varphi }(N^{\ast })={\varphi }_0\varphi (N),$$

$$k^{\ast }(N^{\ast })=k_0{\varphi }_0^{3}k(N),p^{\ast }=(p_0-p_1)p+p_1,{\mathbf{u}}^{\ast }=\frac{k_0{\varphi }_0^3(p_0-p_1)}{\mu L}\mathbf{u},{{\mathbf{u}}_{\mathbf{r}}}^{\ast }=\frac{Q_c}{L}{\mathbf{u}}_{\mathbf{r}},$$

$$D^{\ast }(N^{\ast })=D_{n}D(N),Q_0^{\ast }=\frac{k_0{\varphi }_0^3(p_0-p_1)}{\mu }{Q}_0,{\sigma }^{\ast }=\frac{4\mu \tau Q_c}{L{R}_0{\varphi }_0}\sigma ,t^{\ast }=\frac{t}{\beta },$$

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

$$\frac{\partial N}{\partial t}=F_1(\sigma )N(1-N)-F_2(\sigma )N+\delta \nabla \cdot \left(D(N)\nabla N\right),$$

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

$$\varphi (N)=(1-\rho N),k(N)=\varphi (N{)}^3,\sigma =\frac{\Vert {\mathbf{u}}_{\mathbf{r}}\Vert }{\varphi (N)},$$

$$F_1(\sigma )=1-\left(\frac{1}{2}\right)(\tanh [g_c(\sigma -{\sigma }_c)]+1),F_2(\sigma )=\left(\frac{1}{2}\right)(\tanh [g_c(\sigma -{\sigma }_c)]+1),$$

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

$$p=1\text{at}x=0,0\le y\le 1,\text{and}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).$$

$$q_{ij}^{\ast }=\frac{\pi (R_{ij}^{\ast }{)}^4}{8l\mu }(p_i^{\ast }-p_j^{\ast }),$$

$$\sum _{j=1}^{n^2}{A}_{ij}{q}_{ij}^{\ast }=\{\begin{array}{ll}{f}_{li},& 1\le i\le n,\\ 0,& n<i\le n^2-n,\\ -f_{ir},& n^2-n<i\le n^2,\end{array}$$

$$\sum _{j=1}^{n^2}{A}_{ij}(R_{ij}^{\ast }{)}^4(p_i^{\ast }-p_j^{\ast })=\{\begin{array}{ll}{R}_0^4(p_0-p_i^{\ast }),& 1\le i\le n,\\ 0,& n<i\le n^2-n,\\ -R_0^4(p_i^{\ast }-p_1),& n^2-n<i\le n^2.\end{array}$$

$$Q^{\ast }=\sum _{i=1}^n\frac{\pi R_0^4}{8l\mu }(p_0-p_i^{\ast }).$$

$$q_{ij}^{r\ast }=q_{ij}^{\ast }\frac{Q_l}{Q^{\ast }},$$

$$\frac{d{N}_i^{\ast }}{d{t}^{\ast }}=\beta \left(F_1^l({\sigma }_i^{\ast })N_i^{\ast }\left(1-\frac{N_i^{\ast }}{N_c}\right)-F_2^l({\sigma }_i^{\ast })N_i^{\ast }\right)+\sum _{j=1}^{n^2}{A}_{ij}{D}^{\ast }(N_i^{\ast },N_j^{\ast })(N_j^{\ast }-N_i^{\ast }),i=1\dots n^2,$$

Full text

1 Introduction

Tissue engineering is a rapidly growing field seeking to apply techniques from a variety of disciplines to create tissues and organs. There is increased need for these technologies as patient transplant lists continue to grow, especially due to increased longevity (Beard et al.,, 2013). Despite significant progress in recent years, there is still a need to better understand the basic biological and physical aspects underpinning the engineering of artificially constructed tissues and organs in order to achieve widespread clinical success (Van Blitterswijk and Thomsen,, 2008).

A contemporary strategy for in vitro tissue engineering is to seed cells onto a porous biomaterial scaffold which is placed inside a bioreactor where it is perfused with nutrient-rich culture medium (Glowacki et al.,, 1998; Kim et al.,, 2007; Cimetta et al.,, 2007). There is a large variation in cell types, scaffolds, and bioreactor geometries used to grow artificial tissues. Some cell types are mechano-sensitive and can be induced to proliferate, differentiate, or die in response to the local cellular mechanical environment, e.g. due to the fluid shear stress (Sankar et al.,, 2011; Duan et al.,, 2008; Riha et al.,, 2005; Iskratsch et al.,, 2014; Bakker et al.,, 2004). Some cells produce an extracellular matrix (ECM) that is sufficiently dense to substantially affect the flow of culture medium through individual pores (Hossain et al.,, 2015). Understanding these fluid flow and cell growth interactions, and what roles they have in a tissue engineering experiment, is an important challenge.

Mathematical modelling plays a key role in underpinning experimental design, optimizing operating regimes, and predicting experimental outcomes (O’Dea et al.,, 2012). Mathematical models have also been used in lieu of costly and time-consuming experimental trials. There is a growing literature of models used to understand the interactions between fluid mechanical forces exerted on mechano-sensitive cells, and the effect that cell growth has on the flow via blocking of scaffold pores. The simplest of these are spatially homogeneous models where ODEs are used to model the evolution of each phase, such as in Lemon et al., (2009). However, the spatial structure of a scaffold can have profound effects on the growth of the tissue, due to the transport of cells and nutrients, as well as the variation of mechanical stress (Osborne et al.,, 2010; Lanza et al.,, 2014). Determining the spatial distribution of cells within a porous tissue scaffold, and the nature of their local biomechano-chemical environment, is a key challenge that mathematical modelling can address (Melchels et al.,, 2011; Thevenot et al.,, 2008). Understanding the time-dependent coupling between cell growth and fluid flow is currently an active area of modelling and experimental investigation in tissue engineering and related fields (Geris,, 2013).

To describe the spatial properties of artificial tissue growth, many macroscale continuum models for bioactive porous media have been proposed (O’Dea et al.,, 2012). These approaches are computationally much cheaper than a full simulation of flow within the pores and do not require detailed knowledge of the microstructure. Due to the complexity of growing tissue in vitro, and the wealth of biochemical and biophysical processes involved over a range of spatial and temporal scales, many different theoretical models exist in the literature (German and Madihally,, 2016). These models result in partial differential equations for the dependent variables, such as cell density and fluid velocity, with constitutive assumptions describing the interactions between these variables.

Multiphase models account for multiple phases explicitly (e.g. fluid, scaffold, cells, etc). The governing equations are derived from conservation of mass and momentum for each phase, and interactions between the phases are captured via the specification of appropriate constitutive laws (Pearson et al.,, 2013; O’Dea et al.,, 2010). Alternatively, many models include the effect of cell growth on the fluid flow implicitly, by considering a cell density equation based on conservation of mass, and then prescribing the scaffold porosity to be a function of cell density (Coletti et al.,, 2006; Shakeel et al.,, 2013; Pohlmeyer and Cummings,, 2013; Pohlmeyer et al.,, 2013). We will refer to these as implicit models, in comparison to multiphase models which explicitly account for interactions between phases in terms of physical forces. The flow of fluid through the porous scaffold is then governed by Darcy’s law. In all of these macroscale approaches, constitutive assumptions relate microscale processes, such as cell growth, to macroscale parameters. In both multiphase and implicit models, equations have been derived that account for the influence of shear stress on cell growth. Additionally, homogenization techniques have been proposed to precisely capture the relationship between macroscopic parameters in these models, such as porosity, and growth and flow processes at the pore scale (Shipley et al.,, 2009; Chen et al.,, 2017; O’Dea et al.,, 2015). We briefly give examples of each of these kinds of approaches below.

Pearson et al., (2013); Pearson et al., 2016a ; Pearson et al., 2016b ; Pearson et al., (2015) employed a multiphase approach to model culture medium fluid, scaffold, and cell phases within a hollow fibre bioreactor. This bioreactor design consists of parallel hollow permeable tubes which facilitate nutrient transport to the cells seeded in the extracapillary space surrounding the fibres. The small aspect ratio of the bioreactor was used to simplify these models, and the reduced systems of governing equations were then solved numerically to gain insight into how to stimulate uniform cell growth throughout the scaffold. O’Dea et al., (2010) developed a model of cells, culture medium, and tissue scaffold with the aim of understanding the role of cell-cell and cell-scaffold interactions on tissue growth in a perfusion bioreactor. The focus was on mechanotransduction effects, such as the influence of fluid shear stress on the cell growth. Their analytical and simulation results gave insight into using experimental data to determine the dominant mechanical regulatory mechanisms within a cell population. There is also a very large literature on the use of multiphase models in cancer biology and wound healing, but for brevity we do not review it here; see Byrne and Preziosi, (2003); Byrne et al., (2003) as examples and Roose et al., (2007); Preziosi, (2003); Alarcón, (2009); Powathil et al., (2015); Rieger et al., (2016) for reviews of vascular and avascular tumour modelling. Multiphase models also appear in other areas of bioactive porous material, such as in the development and growth of biofilms (Anguige et al.,, 2006; Klapper and Dockery,, 2010; Fagerlind et al.,, 2012).

We now discuss the class of implicit models. Coletti et al., (2006) proposed a model of nutrient transport, fluid flow, and nutrient-limited cell growth in a tissue-engineering bioreactor which was used to investigate the effect of experimental protocols employed to overcome diffusion-limited transport, and to study adverse conditions that are observed experimentally, such as flow channelling along bioreactor walls. Accounting for both cell density-dependent permeability and mechanotransduction of shear stress by cells, Shakeel et al., (2013) investigated a model of nutrient transport and cell growth in a porous scaffold within a perfusion bioreactor. Cell proliferation and ECM deposition were assumed to decrease the local porosity, and hence permeability of the scaffold, while moderate values of the shear stress were assumed to enhance both cell growth and uptake of a generic nutrient by the cells. Although it was suggested that cell death would occur at higher values of shear stress (corresponding to higher flow rates into the scaffold), this was not explored. Pohlmeyer and Cummings, (2013) and Pohlmeyer et al., (2013) used a similar approach to Shakeel et al., (2013) to model other aspects of cell growth in perfusion bioreactors, including growth-factor driven haptotaxis. Using comparable implicit modelling approaches, Nava et al., (2013) considered a nutrient and shear stress-dependent moving boundary problem to model cell growth in a multi-compartment computational model. Their results included development of computational approaches for coupling sub-models with 2 and 3 dimensional spatial regions. Hossain et al., (2015) considered the coupling of fluid shear stress, cell growth, and the fluid flow in a scaffold-free perfusion bioreactor.

Mathematical homogenization approaches have been proposed to upscale microscale features of the scaffold. Shipley et al., (2009) considered a model for fluid, glucose, and lactate transport in a microstructured porous medium. In this study, the variation of the effective permeability of the scaffold due to cell growth was not considered. O’Dea et al., (2015) considered the advection of a generic diffusible nutrient through a porous medium and derived macroscale equations accounting for microscale accretion of biomass at the pore surface due to nutrient uptake. Homogenisation has also been used to upscale fluid flow and elastic deformation in a fibre-reinforced hydrogel scaffold for cartilage tissue engineering (Chen et al.,, 2017). Outside of tissue engineering, several other homogenisation approaches have been developed to upscale tissue growth processes. Irons et al., (2017) developed a model of biomass accretion in a porous medium aimed at understanding vascular tumour growth. Penta et al., (2014) derived a macroscopic model of a porous linearly elastic medium that accounted for mass exchange at the microscale between the phases due to accretion of the solid phase at the pore surface. Collis et al., (2017) determined effective equations for a poroelastic medium where microscale elastic stresses and deformations were homogenized alongside growth and solute transport. All of the above homogenization results are typically valid for specific parameter regimes, such as when the elastic and growth time scales are strongly separated.

In contrast to the spatially continuous models described above, network models have been used increasingly to describe porous media (Sahimi,, 1994, 1993). These are ‘mesoscopic’ models, whereby the porous medium is discretized as a graph with nodes and edges linking them. Fluid flow can then be described by assigning a pressure to each node, and determining the flow rate between nodes, including any sources or sinks of fluid into the network. Network models of flow through porous media explicitly describing the pore structure have not, to the best of our knowledge, been used in the context of tissue engineering scaffolds. However, network approaches have been used to model flow in geological and geochemical porous media, such as the growth of biofilms in soils or mass accretion in pores due to solid forming reactions, as well as in angiogenesis and its applications to cancer modelling and wound healing (Thullner and Baveye,, 2008; Tsimpanogiannis and Lichtner,, 2012; Anderson and Chaplain,, 1998; McDougall et al.,, 2002). We now review some of these models as motivation for our network modelling approach. We also mention that Mely and Mathiot, (2012); Barbotteau et al., (2003) used percolation theory to model random cell growth in a lattice to determine effective properties of bone tissue engineering scaffolds, although we note that the effect of flow on the growth process was not explicitly modelled, and biomechanical effects such as shear stress were therefore not included.

Biofilms growing in porous media have been modelled using a dynamic (e.g. time-dependent) pore network with coupling between the growth of the biofilm and flow throughout the network. Thullner and Baveye, (2008) proposed a lattice model of a porous medium, where growing biofilms clogged pores. Nutrients were carried by the fluid, affecting the biofilm growth throughout the porous medium. Different constitutive assumptions relating pressure differences to flow rates in pore throats were explored and resulted in global changes in hydraulic conductivity throughout the entire network. Gharasoo et al., (2012) and Rosenzweig et al., (2014) extended this approach to consider the influence of pore heterogeneity, as well as different kinds of reactions, within the network. Tsimpanogiannis and Lichtner, (2012) also considered a reactive transport model utilizing networks in a geochemical setting with solid formation leading to pore blocking, which resulted in changes in the flow throughout the network.

Network modelling has been used extensively in the study of angiogenesis, which is the formation of new blood vessels from existing ones. Alongside vascular network remodelling, angiogenesis plays a major role in cancer biology, wound healing, and developmental biology; see (Otrock et al.,, 2007; Yadav et al.,, 2015; Scianna et al.,, 2013) for contemporary reviews of these fields. The discrete structure of microvessels has motivated a long history of using network models to understand their formation, and associated fluid and mass transport. Anderson and Chaplain, (1998) considered continuum and discrete models of endothelial-cell (tip cell) led vessel sprouting, and specifically derived the (probabilistic) movement of discrete cells via a finite-difference discretization of the corresponding continuum description. Capillaries then form behind these tip cells, growing the vascular network. Extending this approach, McDougall et al., (2002) developed a network model for angiogenesis induced by a nearby tumour, and investigated fluid flow through the resulting network structure. Their results suggested that topological properties of the network, such as the number of anastomoses (loops), were key factors in determining the flow through the resulting network, which impacts the effectiveness of drug delivery via the vasculature. Reviews of some of the earlier theoretical literature utilizing these approaches can be found in Chaplain and Anderson, (2004); Chaplain et al., (2006). More recent network remodelling frameworks have been developed and applied to vascular tumours (Secomb et al.,, 2013; Pries and Secomb,, 2014; Rieger et al.,, 2016; Vilanova et al.,, 2017) whereby vascular remodelling occurs due to lack of flow or tumour angiogenic factors being produced by cells due to hypoxia.

Many recent models of cancer growth and treatment, including many of the network models above, adopt an individual-cell-based modelling approach for some components of the overall tumour model, such as discrete cells. Individual-based models track the local microenvironment around a cell, composed of mechanical stresses and chemical concentrations, in order to capture biophysical interactions between cells and between cells and the surrounding environment (Byrne and Drasdo,, 2009; Fozard et al.,, 2009). These models can either treat the cells as discrete cellular-automata confined to a spatial lattice (Gerlee and Anderson,, 2007), use off-lattice models of cells as points (Drasdo,, 2007), or use cellular-Potts models which specify intracellular and environmental interactions via a Hamiltonian (Glazier et al.,, 2007). Similar discrete modelling approaches have also been employed in vasculogenesis (Perfahl et al.,, 2017). Simulations of many individual cells are then used to gain insight into tissue-level behaviour (Powathil et al.,, 2015).

Comparisons between individual-based and continuum models of growing cell populations have suggested certain parameter regimes where the approaches agree and disagree, which has implications for applicability and computational tractability of each kind of model (Pillay et al.,, 2017; Byrne and Drasdo,, 2009; Osborne et al.,, 2017). Hybdrid discrete and continuum models of vascular tumours have also been studied, where the vasculature itself has a discrete structure, but transport of nutrients, cell migration, and cell/ECM growth throughout the tumour is modelled as spatially continuous (Welter and Rieger,, 2013; Vilanova et al.,, 2017; Welter and Rieger,, 2016; Wu et al.,, 2014; Figueredo et al.,, 2013; Peng et al.,, 2017). These hybdrid models offer the advantage of describing the discrete structure of the vasculature, but can be efficiently simulated computationally (de la Cruz et al.,, 2017). Other multiscale frameworks linking vasculature, cell cycle dynamics, and other effects have also been pursued in the context of infectious diseases and tissue dynamics in general (Cilfone et al.,, 2015; Yu and Bagheri,, 2016), which combine spatially continuous and spatially discrete submodels.

Finally, we also mention analogous approaches to those used in angiogenesis applied to wound healing; these consist of similar continuum, discrete, and hybdrid modelling approaches. McDougall et al., (2012) develops a hybdrid discrete and continuum model of angiogenesis in retinal wound healing, using many of the same ideas of discrete tip cell migration as in Anderson and Chaplain, (1998). Flegg et al., (2015) reviews the current state of reaction-transport models, in particular considering angiogenesis in the context of dermal wounds that heal via replenishment of granulation tissue and ECM. We particularly note the use of discrete models that characterize vasculature, as in the tumour literature discussed above, and continuum models for reaction and transport of cells and nutrients Chaplain and Byrne, (1996); Spill et al., (2015).

Having reviewed some of the network and hybdrid modelling literature, we return to the discussion of a tissue engineering scaffold as an active porous medium. The pore and scaffold length scales in typical experiments with perfusion bioreactors can lead to scaffolds with relatively few pores, and it is unclear that macroscopic spatially continuous models will capture features present in small networks of pores. See Figure 6 of Cox et al., (2015) for an example of a porous scaffold with very few pores, and Loh and Choong, (2013) and Vafai, (2010) for general discussions of pore size, scaffold geometry, and transport phenomena in porous tissue scaffolds. Motivated by concerns with these spatially continuous modelling approaches, we propose a novel lattice model of a bioactive porous medium. We compare this to a typical continuum model, modified from Shakeel et al., (2013), to elucidate the different kinds of behaviour displayed by each modelling paradigm.

We consider the following biological system that captures many of the features encountered when engineering artificial tissues. We model a cell phase, which incorporates the ECM produced by the cells, and a fluid phase modelling the culture medium, which together saturate a two-dimensional porous scaffold. We consider model behaviours on the timescale of cell proliferation, focusing in particular on the role of shear stress on cell proliferation. In Shakeel et al., (2013) and Chapman et al., (2014) it was assumed that moderate levels of shear stress would enhance cell proliferation, but the detrimental effects of large shear stress were not considered. Several studies have shown detrimental effects of high shear on renal epithelial cells (Duan et al.,, 2008) and pancreatic endothelial cells (Sankar et al.,, 2011). The influence of shear stress on cell adhesion and ‘wash-out’ has been studied for human fibroblast cells (Lu et al.,, 2004; Korin et al.,, 2007). We are interested in the effects of high levels of shear stress not explored in previously-mentioned models (O’Dea et al.,, 2010; Shakeel et al.,, 2013; Chapman et al.,, 2014), and assume that cells die if the local shear stress exceeds a certain threshold, either due to specific mechano-transduction mechanisms for apoptosis, or due to detachment induced by the flow which has been investigated in the tissue engineering literature (McCoy et al.,, 2012; Whited and Rylander,, 2014). Under this latter mechanism, our model assumes that cells do not reattach to the scaffold downstream after detachment, which is appropriate for some cell types under high velocity viscous flows (Sircar and Roberts,, 2016). High levels of shear-stress were assumed to induce cell death and detachment in Chapman et al., (2017), where nutrient limitations and other biomechanical stimuli were considered in a hollow-fibre bioreactor. We assume that the culture medium is pumped into the scaffold at a fixed flow rate, and that local increases in cell density decrease the local permeability of the scaffold. For the sake of simplicity, we assume the cells are in a nutrient-rich environment and do not explicitly model a nutrient phase; we instead emphasize qualitative differences between lattice and continuum models with cell growth mediated only by fluid shear stress. Focusing on a single stimulus allows us to elucidate paradigmatic differences between spatially continuous and spatially discrete models in a simple setting.

We note that the lattice model we derive has a similar discrete network structure to those found in the angiogenesis literature, much of which originated in Anderson and Chaplain, (1998), but our model differs in two key ways from the majority of models discussed above. Firstly, the underlying lattice topology is intended to represent the pore space of the tissue scaffold as a static network, rather than evolving a dynamic network. Secondly, we do not track individual cells with random motility, but instead consider nodal cell densities which diffuse via a continuous-time deterministic process. The active medium in our models consists of cells changing the local pore size; cells do not grow or remodel the pore network as is considered in the angiogenesis literature.

In Section 2 we present our continuum and lattice modelling approaches, and discuss the physical and biological parameters. In Section 3 we compare results from these two models, primarily using numerical simulations. Finally in Section 4 we discuss the implications of our study for the successful modelling of porous media in tissue engineering applications.

2 Motivation and Modelling Frameworks

We consider a two dimensional rigid porous medium as a model of the tissue scaffold. For simplicity, we consider a square domain with side length $L$. The culture medium is an incompressible viscous Newtonian fluid, with dynamic viscosity $\mu$. We assume that fluid is pumped at a constant flow rate $Q_{c}$ ($\mathrm{m}^{2}\mathrm{s}^{-1}$) into the left side of the domain, and exits through the right side. We assume no fluid enters or leaves through the horizontal boundaries of the domain. See Figure 1(a) for a visualization of the domain. We assume that cells grow logistically, diffuse within the scaffold, and cannot leave the domain. We assume that cells will stop growing and will die if the shear stress exceeds a threshold value. The timescale of interest is that of cell proliferation. Finally, we neglect advection of the cells by the fluid, as we assume cells which are detached from the scaffold do not reattach downstream.

2.1 Continuum Model of a Bioactive Porous Medium

We describe the square domain as spatially continuous with Cartesian coordinates $\bm{x}^{*}=(x^{*},y^{*})\in[0,L]^{2}$ (we use asterisks throughout to denote dimensional variables). We assume that the pore Reynolds number is small and neglect fluid inertia. We use Darcy’s Law to model flow through the scaffold (Bear,, 1972), so that

[TABLE]

where $\bm{u}^{*}$ is the Darcy velocity of the fluid, $p^{*}$ is the fluid pressure, and $k^{*}(x^{*},y^{*},N^{*})$ is the permeability of the scaffold which depends on the cell density $N^{*}$.

For ease of computation, we prescribe a pressure drop across the domain so that, along with no flux conditions at the horizontal boundaries, we impose

[TABLE]

[TABLE]

where $p_{0}$ is the upstream pressure, $p_{1}$ is the downstream pressure, and $\bm{n}$ is the outward unit normal.

We use the linearity between pressure and the Darcy velocity to rescale the fluid variables to match the prescribed fluid flow rate through the scaffold. The total flow rate through the boundary at $x^{*}=0$ is

[TABLE]

We define the rescaled Darcy velocity to be,

[TABLE]

The cell density $N^{*}$ is governed by the reaction-diffusion equation,

[TABLE]

where $t^{*}$ is time, $N_{c}$ is the maximum cell density, $\beta$ is the cell proliferation rate, $D^{*}$ is a nonlinear cell diffusion coefficient, $\sigma^{*}$ is the fluid shear stress, and the dimensionless functions $F_{1}(\sigma^{*})$ and $F_{2}(\sigma^{*})$ capture the effect of fluid shear stress on growth and death respectively. The terms on the right hand side correspond to logistic growth, cell death, and cell diffusion respectively. For simplicity we assume that $\beta$ also quantifies the death rate due to high shear, and $F_{1}$ and $F_{2}$ are dimensionless functions of the shear stress that vary between [math] and $1$ (see 2.12).

The no flux conditions for the cell density are,

[TABLE]

where $\partial[0,L]^{2}$ denotes the boundaries of the square domain. Finally we impose the following initial condition for cell density,

[TABLE]

where $N_{0}^{*}(x^{*},y^{*})$ is the initial scaffold seeding density.

We now specify how the scaffold permeability $k^{*}$ appearing in Equation (2.1) depends on the cell density $N^{*}$. Following Coletti et al., (2006), we write the porosity as

[TABLE]

where $\phi_{0}$ is the volume fraction of pore space in the scaffold without cells, and $\nu$ is the volume of pore space occupied by an individual cell. We require that $\nu N_{c}\leq 1$ so that under the bounded dynamics of Equation (2.6), $0\leq\hat{\phi}(N^{*})\leq\phi_{0}$ for all $N^{*}\leq N_{c}$.

We follow Shakeel et al., (2013) and relate the cell-dependent porosity to the scaffold permeability via

[TABLE]

where $k_{0}$ is the permeability of the cell-free scaffold.

We require an expression for the shear stress $\sigma^{*}$ in Equation (2.6). We follow Whittaker et al., (2009) and Shakeel et al., (2013) and write,

[TABLE]

as the shear stress experienced by cells at the pore scale in terms of the Darcy velocity, where $R_{0}$ is the typical radius of a pore and $\tau$ is the tortuousity of a typical fluid path. As in Whittaker et al., (2009) and Shakeel et al., (2013), we neglect the influence that cells have on the tortuousity $\tau$ and treat it as a nondimensional constant corresponding to the average ratio of streamline lengths to the straight-line distance between two points. The more tortuous the scaffold, the faster the interstitial fluid has to be to travel longer paths in the same time.

Finally we specify the functions $F_{1}(\sigma^{*})$ and $F_{2}(\sigma^{*})$ to be

[TABLE]

where $g$ and $\sigma_{t}$ are sharpness and threshold parameters. Caricatures of these functions are plotted in Figure 2. These functions model smoothed step-function behaviour and are commonly used in the literature (Coletti et al.,, 2006; Shakeel et al.,, 2013). These functions capture logistic growth for small and moderate shear stress, and cell death at large values of shear stress.

We nondimensionalise as follows

[TABLE]

[TABLE]

[TABLE]

leading to the following nondimensional system of equations,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $\rho=\nu N_{c}$ is the maximum available fraction of the pore space the cells can occupy, $\delta=D_{n}/\beta L^{2}$ is the ratio of proliferation and diffusion timescales, $g_{c}=(4\mu\tau Q_{c})/(LR_{0}\phi_{0})g$ is a sharpness parameter, and $\sigma_{c}=(LR_{0}\phi_{0})/(4\mu\tau Q_{c})\sigma_{t}$ is the threshold parameter.

The nondimensional boundary conditions are

[TABLE]

[TABLE]

[TABLE]

with the initial data

[TABLE]

Equations (2.13)-(2.18) are a modified form of those used by Shakeel et al., (2013), where we have neglected nutrient transport and nonlinear cell diffusion, but included a mechanism for cell death induced by high shear stress.

2.2 Lattice Model

We idealize the scaffold domain of length $L$ as an $n$ by $n$ square lattice representing the connectivity between pores. At each node we specify the local cell density together with fluid pressures which we use to define volumetric flow rates along the edges between nodes. We assume these edges represent pipes of length $l=L/n$. See Figure 1 for a diagram of this idealization.

We prescribe a constitutive law relating the pressure $p_{i}$ at each node $i$ with the volumetric flow rate between the nodes, where $1\leq i\leq n^{2}$. The particular form of this relationship depends on the microstructure of the scaffold under consideration, but for simplicity we assume Poiseuille flow between each node, so that

[TABLE]

where $q_{ij}^{*}$ is the volumetric flow rate ($\mathrm{m}^{3}\mathrm{s}^{-1}$) from node $i$ to node $j$, and $R_{ij}^{*}$ is the effective radius of the pipe between these nodes, which will depend on the cell density at the nearby nodes (see Equation (2.26)).

Due to incompressibility, we can write conservation of mass at each node as

[TABLE]

where $A_{ij}$ is the unweighted undirected adjacency matrix of the lattice representing the connectivity of nodes, and $f_{li}$ and $f_{ir}$ are the volumetric flow rates into the nodes at the left and the right side of the scaffold respectively, which will be used to prescribe the total volumetric flow rate into the scaffold. We have that $A_{ij}=1$ if there is an edge between node $i$ and node $j$, and otherwise $A_{ij}=0$. While we only consider a square lattice here, this formulation is general and any pore network topology can be accounted for by providing a different adjacency structure.

As in Section 2.1 we prescribe a pressure drop across the domain. We fix the pressure upstream of the left boundary as $p_{0}$, and the pressure downstream of the right boundary as $p_{1}$. Then the volumetric flow rate of fluid entering each node at the left is given by $f_{li}=(\pi R_{0}^{4}(p_{0}-p_{i}^{*}))/(8l\mu)$ for $1\leq i\leq n$, and the volumetric flow rate of fluid leaving each node along the right by $f_{ir}=(\pi R_{0}^{4}(p_{i}^{*}-p_{1}))/(8l\mu)$ for $n^{2}-n<i\leq n^{2}$, where $R_{0}$ is the radius of pipes entering and exiting the scaffold.

We combine equations (2.19) and (2.20) together with the expressions for the boundary volumetric flow rates $f_{li}$ and $f_{ir}$, to find that for all nodes $i$,

[TABLE]

Equations (2.21) represent an algebraic system of $n^{2}$ equations for the nodal pressures, $p_{i}^{*}$. As before, the physically relevant boundary condition is a constant volumetric flow rate into the scaffold. We exploit the linearity of the fluid problem above to rescale the fluid variables to match this condition. The total volumetric flow rate into the pipes along the left boundary is

[TABLE]

We rescale the volumetric flow rates $q_{ij}^{*}$ by,

[TABLE]

where $Q_{l}$ is the volumetric flow rate into the scaffold. Note that $Q_{l}$ has dimensions of $\mathrm{m}^{3}\mathrm{s}^{-1}$, whereas for the 2-D spatially continuous model, the flow rate $Q_{c}$ had dimensions of $\mathrm{m}^{2}\mathrm{s}^{-1}$. If we assume that the mean fluid velocity (not the Darcy velocity) into the pore space of the scaffold along the left boundary is the same between the two models, we must have $Q_{c}/(L\phi_{0})=Q_{l}/(n\pi R_{0}^{2})$ where the factor of $\pi R_{0}^{2}$ is the cross-sectional area of each pipe, and $n$ accounts for the number of pipes along the leftmost boundary. Hence $Q_{l}=(n\pi R_{0}^{2}Q_{c})/(L\phi_{0})$, and we use this relationship in the remainder of the paper.

Again assuming that cells grow logistically, die, and diffuse to neighbouring nodes, the evolution equation for the cell density $N_{i}^{*}$ at node $i$ is

[TABLE]

where $D^{*}(N_{i}^{*},N_{j}^{*})$ is the local nonlinear cell diffusion rate to move between neighboring nodes, $\sigma_{i}^{*}$ is a local average of the shear stress defined in Equation (2.29), and the functions $F_{1}^{l}$ and $F_{2}^{l}$ model the effect of shear stress on cell growth and death. Equations (2.24) are analogous to Equation (2.6). In particular, the cell proliferation rate $\beta$ and maximum cell density $N_{c}$ have the same meaning as in Equation (2.6). The cell diffusion rate between nodes can be related to cell diffusion on the length scale of the scaffold $D_{n}$ by $D_{l}=D_{n}/l^{2}=n^{2}D_{n}/L^{2}$. We specify an initial condition at each node as

[TABLE]

We model the effect of cell growth on the fluid flow by taking the effective radius of a pipe to depend on the nearby cell densities as

[TABLE]

where again $\nu$ represents cell volume. Equation (2.26) models pipes with radii that linearly decrease as the cell density increases at nearby nodes. Note that we must have $\nu N_{c}\leq 1$ for this radius to be non-negative for all feasible cell densities. This relationship is analogous to the porosity relationship given in Equation (2.9).

We compute the shear stress in each pipe under the assumption of Poiseuille flow made in Equation (2.19). The local velocity profile in the pipe $R_{ij}^{*}$ is

[TABLE]

where $0\leq r^{*}\leq R_{ij}^{*}$ is the radial coordinate. The magnitude of the shear stress at the wall of each pipe is then

[TABLE]

As cell densities are defined at the nodes, we must relate the shear stress in each idealized pipe to an average shear stress $\sigma_{i}^{*}$ at each node. We note that adding the shear stresses from the four neighboring nodes accounts for the total volumetric flow rate twice, so we sum the shear stress in each pipe connected to node $i$ and divide by $2$ to obtain

[TABLE]

We model the influence of this averaged shear stress on cell proliferation and death as before by

[TABLE]

These functions have the same form as those in Equations (2.12); see Figure 2 for a visualization.

For each $1\leq i,j\leq n^{2}$ we nondimensionalize by taking

[TABLE]

[TABLE]

from which we obtain

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $i=1\dots n^{2}$, $\delta=(D_{n})/(L^{2}\beta)$ is the ratio of proliferation and diffusion timescales, $\rho=\nu N_{c}$ is the maximum available fraction of pore radius cells can occupy, $g_{l}=(2\mu Q_{c}g)/(LR_{0}\phi_{0})$ is a sharpness parameter, and $\sigma_{l}=(LR_{0}\phi_{0}\sigma_{t})/(2\mu Q_{c})$ is the shear stress threshold parameter. Note that $\rho$ and $\delta$ are identical parameters between the lattice and continuum models. The parameters $g_{l}$ and $\sigma_{l}$ are analogous to $g_{c}$ and $\sigma_{c}$ but may have different values for a particular experimental system. The number of pores, $n$, does not explicitly appear in the cell density equation (except to scale diffusion) or the constitutive law for its relationship to pipe radii as we have accounted for issues of lengthscale in thinking of $N_{i}$ as the nondimensional density of cells, rather than a cell number. We also have the initial data

[TABLE]

In addition to fundamental differences between a spatial continuum and a discrete lattice, there are two significant constitutive differences between our models. The exponents in the relationships between the Darcy velocity, $\bm{u_{r}}$, or volumetric flow rates, $q_{ij}$, which play analogous roles, and the cell density are not the same between the two models. To see this, we write Equation (2.13a) using the permeability $k$ and porosity $\phi$ in Equations (2.13f,g) as

[TABLE]

and similarly rewrite (2.31a) using the effective pipe radii from Equation (2.31f) to find,

[TABLE]

These relationships show that for a given local pressure drop, the local Darcy velocity in the PDE is larger in regions of high cell density than the corresponding volumetric flow rate in the lattice.

There are also differences in the nondimensional expressions for the shear stress in each model (see Equation (2.13h) and Equation (2.31g)). To visualize how these constitutive differences lead to quantitatively different predictions, consider a uniform cell density in both models. Let $N(x,y,t)=\hat{N}(t)$ for all $(x,y)$ in the domain and $N_{i}(t)=\hat{N}(t)$ for $1\leq i\leq n^{2}$. For a uniform cell density, the fluid flow is uniform so $\|\bm{u_{r}}\|=1$ and $q_{ij}^{r}=1/n$ for all $i,j$. This can be seen from solving Equations (2.13a,b) and (2.13d-h) for the PDE, or Equations (2.31a,b), and (2.31d-g) for the lattice. Note that in the lattice, we have nondimensionalized such that the total volumetric flow rate into the scaffold from the left boundary is $1$, leading to the flow in each horizontal pipe being $1/n$ as there are $n$ pipes along the boundary. We can then compute the (spatially uniform) shear stress for the PDE model as,

[TABLE]

and for the lattice model as,

[TABLE]

We plot these shear stresses as functions of $\hat{N}$ in the interval $[0,1]$ in Figure 3. For large cell densities, the shear stress experienced by the cells in the continuum and lattice models differs by two orders of magnitude. These differences between Equations (2.36) and (2.37) prevent quantitative agreement between our models; we emphasize qualitative differences between the model predictions by exploring ranges of parameters in the following discussion.

Various forms of nonlinear cell diffusion have been used to model the movement of cells within porous media. For simplicity, we primarily consider only linear cell diffusion ($D(N)=D(N_{i},N_{j})=1$), so that both models have linear cell diffusion at a rate $\delta$. In Section 3.5 we demonstrate that the linear diffusion regime still accurately captures the main features we are trying to illustrate.

2.3 Model Paremeters

The continuum and lattice models contain a number of parameters that can vary significantly between tissue engineering experiments. As we are primarily interested in informing model selection, rather than constructing quantitatively accurate models, we focus on the orders of magnitude over which these parameters vary. Table 1 lists typical parameter values found in the literature. Using these values we estimate the ranges of the two nondimensional parameters, $\rho$ and $\delta$, that appear in both models, and estimate the number of pores per side of the square lattice $n$.

The shear stress threshold $\sigma_{t}$ and the sharpness parameter $g$ are not readily available in the literature. We treat the nondimensional thresholds $\sigma_{l}$ and $\sigma_{c}$ as model parameters and demonstrate model behaviours as they vary. Motivated by existing theoretical models in which the function $F_{1}$ approximated a step function (Shakeel et al.,, 2013), we set $g_{l}=g_{c}=60$. In our simulations we fix $\rho=0.9$ so that cell growth significantly affects the effective permeability of the scaffold.

The nondimensional parameters $\rho$ and $n$ can readily be determined from the properties of a specific porous scaffold microstructure and cell type. We can use the values in Table 1 to estimate that the number of pores in a typical scaffold ranges from approximately $10^{2}$ to $10^{6}$. Similarly, given the values of the cell volume $\nu$ and pore radii $R_{0}$, we compute that the number of cells required to fill a pore can vary in the range of $10^{2}-10^{4}$. We now simulate our models in order to explore behaviours as these nondimensional parameters vary.

2.4 Numerical Simulations

The lattice size $n$ and diffusion to proliferation rate $\delta$ are varied to demonstrate the range of qualitative behaviours displayed by these models. Specifically we take $n=25,50,75$, and $100$ for the lattice simulations, and $\delta=10^{-4},10^{-3},10^{-2}$, and $10^{-1}$ for both models. We take values of both lattice and PDE thresholds to be $\sigma_{c}=\sigma_{l}=2.5,5,7.5,10,100$, and $1000$.

The continuum model given by Equations (2.13)-(2.18) was simulated using the finite element solver Comsol with 24,912 triangular elements. Time and space refinements were carried out to ensure that the numerical approach converged. Additionally, the results found via Comsol were consistent with those found from a finite-difference scheme implemented to solve the same model. The lattice model given by Equations (2.31)-(2.33) was solved using an explicit adaptive Runge-Kutta method in Matlab. To ensure accuracy of our simulations with respect to the bifurcation behaviour discussed in Section 3.3, we constrain the maximum Runge-Kutta time step to be $10^{-3}$ (Christodoulou,, 2008). Specific simulations were also undertaken with a fixed time step Runge-Kutta scheme with refinements in the size of the time step to ensure convergence.

We consider our initial condition to be a perturbed uniform cell density. On a lattice with $n=100$ nodes per side, we set $N_{i}(0)=0.1+h\xi_{i}$ where $\xi_{i}\sim\mathcal{N}(0,1)$ is a normally distributed noise term, and we use three values of the variance, $h=10^{-2},10^{-3},$ and $10^{-4}$. Each realization of this perturbed uniform state is then interpolated on the smaller lattices, or onto the triangular elements of the continuum model, so that all simulations have approximately consistent initial conditions. We compute mean cell densities as,

[TABLE]

for the lattice and

[TABLE]

for the PDE (computed numerically using the triangular elements from the finite element scheme). For each set of parameters simulated, and each value of the variance $h$, we repeat the simulation 100 times using different random seeds in order to investigate how robust each result is to different initial cell densities.

3 Results and Discussion

In Section 3.1 we discuss how the cell density changes from uniform logistic growth to death in local regions of the scaffold, leading to spatial heterogeneity in cell density which is observed over long timescales. In Section 3.2 we discuss the onset time of this spatial patterning. In Section 3.3 we discuss oscillations in cell density that are observed only in simulations of the lattice model. In Section 3.4, we compare the model predictions for the final cell densities. Finally, in Section 3.5, we compare the results of a linear cell diffusion model with those of a model with nonlinear cell diffusion.

3.1 Spatial Heterogeneity

Over long timescales (nondimensional times of $t>10$) all of our simulations exhibit significant spatial heterogeneity in cell density and shear stress. In some regions of the scaffold the cell density is close to the nondimensional carrying capacity, and there is reduced fluid flow and low shear stress in these areas. There are also regions of high fluid flow and associated shear, and low cell density. The spatial structure of the cell density distributions depends heavily on the parameter $\delta$.

For large values of the parameter ($\delta\geq 10^{-3}$), a single region or channel develops with low cell density, and through which the majority of fluid passes. In the remainder of the domain, there is little fluid flow and the cell density is high; see Figure 4. The location of the large fluid channel, corresponding to regions of low cell density, depends on the initial data and parameter values, but we always observe a single channel for $\delta\geq 10^{-3}$. For brevity we plot only the cell density distributions in the scaffold, as the associated shear stress can be inferred from these plots (the fluid flow is always low in regions of high cell density and vice versa).

For $\delta=10^{-4}$, we find significantly more spatial structure in the final cell density distributions. Figure 5 shows plots of the cell density at various times for lattice sizes $n=25,50$ and $100$, with threshold parameter $\sigma_{l}=2.5$ (for one realization of the initial cell density). After a period of growth, the initial uniform state rapidly breaks up into several aggregates of high cell density (e.g. $N_{i}>0.8$, typically close to $1$) surrounded by regions of low cell density (e.g. $N_{i}<0.2$, typically close to [math]) surrounded by regions of low cell density, which then evolve slowly due to diffusion and shear stress mediated cell growth and death. For the smallest lattice size ($n=25$), this process stabilizes quickly leading to a steady state cell density within the first $t=10$ time units (compare the cell density distributions for $n=25$ and $t=10,15,$ and $30$ in Figure 5). For $n=100$, many interior clusters slowly shrink as the cell aggregates along the top and bottom boundaries of the scaffold continue to grow over a longer timescale. Figure 6 shows cell density distribution plots for $\sigma_{l}=7.5$. Note that for each $n$ and $t$, the total area of high cell density is larger in Figure 6 than in 5. For both $\sigma_{l}=2.5$ and $\sigma_{l}=7.5$, the number of clusters of high cell density is initially greater in the $n=100$ case after the onset of non-uniform growth, but there is a coarsening over time which reduces the number of distinct fluid channels and disconnected regions of high cell density. In comparison, the number of disconnected regions of high cell density is larger in the smaller lattices by the end of the simulation time in both Figures 5 and 6 (compare the $t=30$ plots for $n=25$ and $n=100$ in both Figures). This difference in the number of disconnected regions of high cell density between small and large lattices is robust across multiple realizations.

Note that in Figure 5 the coarsening for the $n=100$ lattice involves small clusters of high cell density disappearing, while for Figure 6 it is primarily channels of low cell density that disappear over time. This difference in long-time behaviour can be understood via the threshold parameter $\sigma_{l}$. When $\sigma_{l}=2.5$, cells die at relatively low values of shear stress and the overall cell density is low so cells recede away from regions of high fluid flow. For $\sigma_{l}=7.5$ cells can withstand much more shear and grow over a much larger proportion of the scaffold, filling in small channels of fluid flow over time.

These results demonstrate a diverse range of behaviour for the lattice models. Figure 7 illustrates predictions from the PDE model for the cell density distribution over the spatial domain at various times and for four values of the threshold parameter $\sigma_{c}$ for one realization of the initial cell density. This model also exhibits transient dynamics where clusters of high cell density form, and then coarsen into larger regions of high cell density as in the $n=100$ lattice shown in Figures 5 and 6. Overall there are more regions of high cell density in the PDE model compared to the lattice simulations, which is consistent with the shear stress predictions given in Figure 3 (for the same cell density, this plot shows lower predicted shear stress values for the PDE model compared to the lattice model).

In Figure 8, we plot time series of mean cell densities, $\hat{N}(t)$, for lattice simulations with $n=25,50,$ and $100$, as well as for the PDE for the same realization of the initial cell densities used in Figures 5-7. The $n=25$ lattice quickly reaches equilibrium for each $\sigma_{l}$, whereas the larger $n=100$ simulations show growth and death processes over a longer timescale, and in some cases these have not reached an equilibrium value at the end of the simulations at $t=30$ (compare the long time behaviour in Figures 8(a) and 8(c)). We truncate the numerical experiments here partly because temporal variations in mean cell density beyond this point for all simulations were small, and partly because this period of time would exceed most in vitro tissue engineering experiments as $t=30$ would typically be on the order of months. All simulations shown in Figure 8 grow logistically at the same rate until some region of the scaffold reaches the shear stress threshold $\sigma_{l}$ or $\sigma_{c}$ leading to non-uniform growth, and eventually to differences in final mean cell densities.

We observe less spatial structure in the PDE simulations than in the lattice simulations (compare Figures 5 and 6 with Figure 7). We quantify spatial variation in cell density with a heterogeneity function. For simplicity we choose

[TABLE]

which is the total absolute difference in cell density between adjacent nodes throughout the lattice. The function satisfies $E\geq 0$ with $E=0$ only for a homogeneous cell density distribution. The PDE solutions were interpolated onto a square grid with a size of $100$ by $100$, so the function defined by (3.1) can be seen as a discretization of the functional $\int_{0}^{1}\int_{0}^{1}\left|\partial_{x}N(x,y)\right|+\left|\partial_{y}N(x,y)\right|dxdx$ which is a measure of anisotropic total variation (van Gennip and Bertozzi,, 2012). This heterogeneity function allows for a quantitative visualization of the coarsening behaviour described earlier that was qualitatively observed in Figures 5-6.

In Figure 9 we plot values of this function for the simulated cell densities at $t=30$ averaged over $300$ realizations of the initial cell densities. For the majority of parameter combinations simulated, the value of the heterogeneity given by (3.1) was lower for the PDE simulations than for any of the lattice simulations with the same parameters. Larger values of diffusion show significantly smaller values of spatial heterogeneity, and for most combinations of the parameters the larger lattice of size $n=100$ has a lower value of $E$ than the other two lattice sizes. Figure 10 shows a time series of this function for one choice of parameters for the lattice simulations and for a corresponding PDE simulation. For $n=25,50$, and for the PDE, the difference between the final heterogeneity and its maximal value is small, while for $n=100$ the maximum value of $E$ occurs around $t\approx 4$, and then slowly falls over time to a value that is significantly smaller (compare this to Figure 6 for $n=100$ and $t\geq 5$).

There is less spatial heterogeneity in the final cell densities of the PDE simulations than in the corresponding lattice model for all values of $n$ and all values of $\sigma_{l}$ and $\sigma_{c}$ except for $\sigma_{l}=\sigma_{c}=2.5$, where we note that for the cell density distributions in Figure 5 the cells have aggregated to the sides for $n=100$ and $t=30$ and hence created a single large channel in the center of the scaffold, whereas several small channels remain in Figure 7. We note that in Figure 10, there is a specific time at which spatial heterogeneity emerges (e.g. $E$ is no longer [math]), which corresponds to the onset of cell death in some regions of the scaffold.

3.2 Onset of Non-Uniform Patterning Due to Shear-Induced Cell Death

As the initial data are approximately uniform, the evolution of the cell density can be determined analytically up to the time that cell death occurs. For $\sigma\lesssim\sigma_{c}$, the cell growth is purely logistic. This bound is sharp in the limit of $g\to\infty$ where $F_{1}$ and $F_{2}$ become step-functions, but it is approximately obeyed for $g=60$ so that for shear stress values below $\sigma_{c}$, cells will grow logistically. Hence, for a constant initial condition $N_{0}\in[0,1]$, we have

[TABLE]

where $t_{s}$ is the time that it takes for the shear stress to become approximately equal to the threshold parameter, and hence for the cells to no longer undergo uniform logistic growth everywhere.

From Equation (2.36) we see that $t_{s}$ may be defined implicitly by

[TABLE]

Substituting (3.2) into (3.3) gives

[TABLE]

Analogously for the lattice model, for $\sigma_{i}\lesssim\sigma_{l}$, the cell growth is purely logistic at each node. So for a constant initial condition $N_{i0}=N_{0}\in[0,1]$ for each $i$, we have

[TABLE]

$t_{s}$ again being the time for the shear stress at some node to become comparable to the shear threshold. From Equation (2.37) we approximate $t_{s}$ by

[TABLE]

We substitute (3.5) into (3.6) to find

[TABLE]

We compare Equations (3.4) and (3.7) to the numerically computed values of $t_{s}$ which correspond to $t_{s}\equiv\min_{t}\left(\max_{x}\left(\sigma(\bm{x},t)\right)>\sigma_{c}\right)$ for the PDE simulations, and $t_{s}=\min_{t}\left(\max_{i}(\sigma_{i}(t))>\sigma_{l}\right)$ for the lattice simulations. In Figure 11 we plot values of $t_{s}$ for all lattice and PDE simulations with varying $\delta$, $\sigma_{c}$, and $\sigma_{l}$. We see excellent agreement between the analytical and numerical predictions for both the lattice and the PDE. We see that $t_{s}$ does not depend on the lattice size $n$ or the parameter $\delta$ except for the PDE simulation at $\sigma_{c}=10$. We note that $\sigma_{c}=10$ is a limiting case as Equation (3.4) indicates that $t_{s}\to\infty$ as $\sigma_{c}\to 10$ if the initial cell density distribution was exactly uniform.

The time $t_{s}$ can be seen in Figure 8 where the form of the time series changes from a logistic curve, and for all values of $\sigma_{l}$ and $\sigma_{c}$ the mean cell density begins decreasing at this time. In most cases this period of decreasing mean cell density is short relative to the simulation timescale. The time $t_{s}$ can also be seen in Figure 10 as the point where the heterogeneity increases sharply from $E=0$, due to the onset of non-uniform growth.

3.3 Lattice Oscillations

A particularly interesting behaviour we observe in some lattice simulations are oscillations in cell density. Such oscillations are not found the PDE model. Figure 12 shows a time series plot of the cell density at every lattice node for $n=10$ and $n=25$ in (a)-(b), and corresponding plots of nodal values of the shear stress in (c)-(d). After a period of transient behaviour, the nodal cell densities and nodal values of the shear stress oscillate in phase. Nodes that are growing ($\sigma_{i}<\sigma_{l}$) are green solid lines, and nodes that are dying ($\sigma_{i}>\sigma_{l}$) are red dashed lines. The spatial mean cell density and shear stress are plotted in blue, which is also oscillating but with a small amplitude. Thus the overall effect of the lattice oscillations on the mean cell density is small. The particular time-series shown in Figure 12 will have different transient behaviour for different realizations of the initial data, but due to the simple spatial behaviour at large diffusion parameter (e.g. Figure 4), the amplitude and frequency of the oscillations for long times is the same for each realization we simulated.

In regions of the scaffold where the local shear stress exceeds the threshold $\sigma_{l}$, cells are (exponentially) dying due to high shear stress. These nodes all have relatively low cell density and hence more fluid passes through them, maintaining this high value of shear stress. In other regions of the scaffold with low shear stress, cells grow logistically. Diffusion acts to move cells from regions of high density to low density, and so cells move from regions with low shear stress to regions of high shear stress. In certain parameter regimes, this process finds an equilibrium value where these regions separate, as in Figures 5 and 6, where between regions of growth and death there are some nodes of intermediate cell density. However, for certain values of $n$, $\sigma_{l}$, and $\delta$, the movement and growth of cells gives rise to a different behaviour from this equilibrium as the growth of cells in regions of high cell density substantially affects fluid flow and shear stress throughout the entire scaffold. As seen in Figure 12, diffusion from growing regions can cause regions that are dying to increase in total cell density, and similarly, regions that are growing logistically can decrease in cell density due to diffusion.

We plot bifurcation diagrams of the maximal nodal oscillation amplitudes in Figures 13(a)-13(b) for several lattice sizes and two values of the threshold parameter, over a range of the parameter $\delta$. After enough time has passed to ensure we are no longer observing transient dynamics, we compute the maximum nodal amplitude as $N_{osc}=\max_{i}\{\max_{t}\{N_{i}(t)\}-\min_{t}\{N_{i}(t)\}\}$ to capture the largest oscillation. We plot the frequency $\omega_{N}$ of these nodal oscillations computed with the Fast Fourier Transform in Figures 13(c)-13(d). We note that the frequency is the same for every node in the lattice.

These bifurcation diagrams show that for each $n$ and $\sigma_{l}$, the one-dimensional parameter space over $\delta$ is composed of disconnected regions where oscillations are permitted (see behaviour illustrated in Figure 12), and regions where cell densities tend to steady states (see Figure 8). Further simulations within each of these regions (using different realizations of the perturbations of the initial data, for instance) have consistent behaviour, and we conjecture that that these behaviours are generic within each region in parameter space demarcated by the bifurcation between steady state behaviour and oscillations. In Figure 13 we see that the number of disconnected regions in $\delta$-space with oscillatory behaviour increases as $n$ increases, and the maximal amplitudes decrease with increasing $n$. We note that for $n=100$, there are many disconnected bifurcation regions for very small $\delta$, but these all have correspondingly small magnitudes in oscillation and so they are not visible in Figure 13. For larger values of $n$, the oscillation amplitudes become so small that they are comparable to numerical truncation and rounding errors, and we conjecture that these oscillations play no role for larger lattices. There is a trend of increasing oscillation frequency and decreasing nodal amplitude for larger values of $\delta$ for a given lattice size. The larger lattices have significantly smaller variations in mean cell density of the scaffold due to these oscillations. This is due to only a small subset of the nodes switching between growth and death behaviours during an oscillation. These oscillations do not occur in any of our PDE simulations. For large $n$, our simulations show that the amplitudes of oscillation in mean cell density decreases, as shown by the blue lines in Figure 12, and that the proportion of nodes switching between growth and death behaviours also decreases. Similarly, nodal shear stress oscillations also decrease in amplitude as the lattice size increases. In the other direction, for lattices of size $n=5$ and smaller, no oscillations were observed. We refer to (Krause et al.,, 2017) for a discussion of Hopf bifurcations in a similar lattice system, where cells die due to high fluid pressure rather than shear stress. There we argue that the oscillations observed in this model are due to a combination of nonlocal effects from the quasi-static fluid equations, coupled with symmetry-breaking bifurcations in the lattice models.

3.4 Mean Cell Density Predictions

We can broadly understand the global behaviour of the numerical simulations by comparing their final mean cell densities computed via Equations (2.38) and (2.39). Figure 14 is a plot of the spatial mean of the cell density after $30$ nondimensional units of time. For any given shear stress threshold and $\delta$, the PDE model has a higher mean cell density. This is expected due to the differences between the constitutive assumptions for shear stress in each model (see Figure 3).

The value of the parameter $\delta$ has a non-monotonic effect for some lattice simulations. In particular, for $\sigma_{c},\sigma_{l}\leq 10$, the maximum mean cell density for the $n=50$ lattice occurs for the intermediate value of the diffusion parameter $\delta=10^{-3}$. For the large value of $\delta=10^{-1}$ in each plot, the lattice size is insignificant in determining the mean cell density. Similarly, for the large shear stress thresholds $\sigma_{l}=100$, and $1000$ shown in the bottom two plots, cell death occurs only in a very small region of the domain, so the lattice size is almost inconsequential in determining the final mean cell density. For $\sigma_{l}\leq 10$ and $\delta\leq 10^{-2}$, however, there are differences between the final mean cell density for different lattice sizes. This suggests that for applications where the ratio of proliferation to diffusion timescales is small, and where cells are very sensitive to high shear stress, the topology of the underlying pore network plays a role in determining the final mean cell density.

3.5 Nonlinear Cell Diffusion

Here we consider the effect of nonlinear cell diffusion. We follow Shakeel et al., (2013) and use the following (nondimensionalized) functional forms for the diffusion coefficient,

[TABLE]

for the continuum and lattice models respectively. For values of cell density near the carrying capacity (dimensionless value of 1), diffusion given by (3.8) approximates the linear case with the same nondimensional parameter $\delta$. For small values of cell density, the effective diffusion of cells throughout the porous medium is smaller than the corresponding linear value ($\delta$) for all $n$. The dimensionless parameter $\gamma$ determines the strength of this nonlinear effect; Shakeel et al., (2013) used $\gamma=2$ to ensure a sufficiently nonlinear behaviour while remaining within a parameter regime where the speed of a proliferating front of cells could be computed. Here we carry out several simulations for varying values of $\gamma$ in the interval $(0,5]$. For small nonlinearity ($\gamma\lesssim 10^{-2}$), nonlinear cell diffusion simulations have identical mean cell density and oscillatory behaviours as linear cell diffusion simulations. For moderate nonlinearity ($10^{-2}\lesssim\gamma\lesssim 3$), qualitative features such as cell density aggregation and oscillations are retained, but quantitative differences in cell density distributions and oscillation frequency appear. For larger values of $\gamma$, we no longer observe oscillations in cell density.

We first consider the spatial cell density distributions with nonlinear cell diffusion. We take $\gamma=2$ in Equations (3.8) and repeat the simulations (using the same realization of the initial data) from Figure 6, and display these in Figure 15. While the steady state cell density distributions have changed, the organization of cells into high density aggregates occurs in both the linear and nonlinear cell diffusion simulations. Additional simulations (with different realizations of the initial cell density, and different values of $\sigma_{l}$) show qualitatively similar coarsening behaviour in the $n=100$ simulations with nonlinear cell diffusion compared to the linear case (see Figures 5-6). Similarly, $n=25$ simulations including nonlinear cell diffusion quickly settle into a steady state as in the linear case. Quantitatively, the change in the spatial mean cell density for $\gamma\leq 1$ is less than 5% in each case shown in Figure 14. Larger values of the nonlinearity parameter $\gamma$ leads to larger changes to the spatial mean cell densities reported in Figure 14 for linear cell diffusion, but we leave further quantification of this to future work.

Next we consider the effect of nonlinear cell diffusion on the oscillations in cell density found in Section 3.3. For small values of the nonlinearity, e.g. $\gamma\leq 10^{-2}$, the oscillatory behaviour of solutions is quantitatively the same up to a slight phase shift due to differences in transient dynamics. We plot two example time series in Figure 16 that use the same parameter sets as those in Figure 12, but now including nonlinear cell diffusion with $\gamma=0.1$ and $1$. We see that the period of the oscillation has approximately doubled in both cases compared with the linear diffusion simulations shown in Figure 12. For increasing values of $\gamma$, many of the regions (in $\delta$) of oscillatory behaviour shown in Figure 13 become smaller, and completely disappear for $\gamma\geq 3$. We conclude that large values of nonlinear diffusion can be stabilizing in that oscillatory regions of parameter space shrink in size, and all oscillatory behaviour disappears for large enough values of $\gamma$. As before, we do not observe any oscillatory behaviour in the continuum model with nonlinear cell diffusion.

4 Conclusions

In this paper we have explored complementary lattice and continuum models for a bioactive porous tissue scaffold, and demonstrated several important differences in model behaviours. The lattice simulations show that considering finite pore networks when modelling cell growth within a porous scaffold will lead to qualitative differences compared with PDE models. Varying degrees of spatial heterogeneity in the cell density distribution are displayed in the lattice model, whereas for most of the parameter space the PDE solutions exhibit simpler spatial structures, such as one or two horizontal channels. The significant variation between different lattice sizes indicates that considering the spatial network explicitly has a nontrivial impact on pattern formation. The smoothing behaviour observed in both the PDE and large $n$ lattice models is reminiscent of models of coarsening in condensed matter physics, where initial spatial irregularities evolve into larger spatial structures. That this coarsening occurs at different timescales and to different degrees due to the size of the lattice is an unexpected result. We conjecture that, modulo constitutive differences, a PDE and a very large lattice (e.g. $n>100$) should have comparable behaviour, but for smaller pore networks the finite structure of the lattice becomes important to the overall growth process within the scaffold.

We observe small amplitude oscillations in cell density exclusively in the lattice model. While the effects of these oscillations on the total mean cell density of the scaffold are small, and would likely not be physically observed given experimental noise, their presence is an interesting effect of finite lattice size. Additionally, the larger variations in nodal values of shear stress (see Figure 12(c)-(d)) has important implications for shear-stress mediated tissue growth. In particular, shear-stress can play a role in both cell proliferation and cell differentiation (Iskratsch et al.,, 2014). These spatiotemporal oscillations in shear stress could potentially lead to the formation of heterogeneity in the cell populations within the scaffold, and hence influence the integrity and uniformity of artificial tissue (Yin et al.,, 2016).

Nonlinear interactions could amplify these oscillations. For instance, some tissue scaffolds are cyclically loaded to facilitate cell growth via mechano-transduction (Angele et al.,, 2004; Pohlmeyer and Cummings,, 2013; Nessler et al.,, 2016). In certain parameter regimes, this oscillatory forcing could interact in non-trivial ways with these lattice induced oscillations, for instance inducing resonance in the cell density oscillations. This is especially true given that small amplitude oscillations in cell density give rise to larger shear stress variations, which suggests feedback between cell proliferation and shear stress. We leave the investigation of such possibilities to future work.

Some strategies for quantifying experimental cell distributions have been proposed which could validate the steady-state distributions of cell density we demonstrate in Figures 5-7 (Thevenot et al.,, 2008), though we note that these approaches involve sacrificing tissue scaffolds after cells have stopped proliferating. Observing oscillations in cell density experiments would be difficult given the small oscillation amplitude, and the difficulty in producing temporal datasets from tissue scaffolds due, e.g., to requiring the scaffolds to be sacrificed in order to ascertain cell density. Models can also be validated by comparing other key outputs, such as the mean cell densities given in Figure 14, against experimentally accessible data (O’Dea et al.,, 2012). While validating these models is difficult, once validated such theoretical models provide predictions for the spatiotemporal evolution of important quantities (e.g. shear stress and cell density) throughout the scaffold, without the need for large numbers of costly and time-consuming experiments.

While we observed quantitative differences between these paradigms (see Figure 14), we attribute some of these differences to the constitutive assumptions as shown in Equations (2.34)-(2.37), and these can vary between particular applications of these models. We have explored simulations using different functional forms of $\phi$, $k$, and $R_{ij}$, as well as different values of $\rho$. These preliminary simulations were consistent with the qualitative differences reported here between lattice and continuum models. The constitutive choices we have made in this paper do give results that are consistent with Nava et al., (2013) in that our 2-D fluid model in the PDE over-predicts cell density compared with the 3-D fluid (‘tube’) model of the lattice. We also note that despite the simplicity of our models, the qualitative differences between lattice and continuum formulations exist even when we introduce nonlinear diffusion into the model. We conjecture that these differences are inherent to the forms of the models used here suggesting that discrete spatial geometry should be carefully incorporated into macroscopic tissue engineering models. This is in line with the study of discrete and continuum cell monolayers (Byrne and Drasdo,, 2009), where insights from individual and lattice cellular-automata models were necessary to inform the continuum formulation.

These differences in the type of behaviours displayed by continuum and lattice models exist even though the network geometry studied here, a square grid lattice, is simple. More complicated hexagonal and face centred cubic lattices (Welter and Rieger,, 2013), off-lattice models (Vilanova et al.,, 2017), and other complex microscale models have been used in hybdrid discrete-continuum settings. In addition to a more complex network geometry, there are many components that could be included in this lattice model, such as a model of nutrient transport, or a more detailed microscale relationship between cell density and fluid flow. Nevertheless, these simple models give insight into the kinds of differences one can expect from each modelling paradigm.

For experiments involving small pore sizes or large constructs, the pore network is relatively dense and we expect continuum and lattice models to be comparable in terms of predicting cell proliferation and global fluid properties, such as scaffold permeability. For small scaffolds or scaffolds with large pores, however, we expect lattice models to significantly differ from continuum models in their predictions. The variations between the lattices of different size implies that a good model of the finite network geometry of a scaffold is important in understanding the tissue growth process for such scaffolds. Given a realistic description of the evolution of cells and fluid flow at the pore scale, our approach gives an alternative to continuum approaches such as mathematical homogenization in upscaling these microscale processes. Numerical solution of the lattice model is no more difficult than an equivalent discretization of a continuum model, and so the only practical disadvantage of using a lattice model would be the requirement of specifying the network topology, and determining appropriate fluid and cell properties at the scale of nodes in the network. We believe that pursuing these kinds of models can lead to novel insights in understanding the growth of artificial tissue, and eventually in developing clinically successful technologies.

5 Acknowledgements

D. Beliaev was partially funded by Engineering & Physical Sciences Research Council (EPSRC) Fellowship ref. EP/M002896/1. We gratefully acknowledge the anonymous referees for providing critical analysis that significantly improved an earlier version of this manuscript.