Curvature- and fluid-stress-driven tissue growth in a tissue-engineering scaffold pore

Cell proliferation within a fluid-filled porous tissue-engineering scaffold depends on a sensitive choice of pore geometry and flow rates: regions of high curvature encourage cell proliferation, while a critical flow rate is required to promote growth for certain cell types. When the flow rate is too slow, the nutrient supply is limited; when it is too fast, cells may be damaged by the high fluid shear stress. As a result, determining appropriate tissue-engineering-construct geometries and operating regimes poses a significant challenge that cannot be addressed by experimentation alone. In this paper, we present a mathematical theory for the fluid flow within a pore of a tissue-engineering scaffold, which is coupled to the growth of cells on the pore walls. We exploit the slenderness of a pore that is typical in such a scenario, to derive a reduced model that enables a comprehensive analysis of the system to be performed. We derive analytical solutions in a particular case of a nearly piecewise constant growth law and compare these with numerical solutions of the reduced model. Qualitative comparisons of tissue morphologies predicted by our model, with those observed experimentally, are also made. We demonstrate how the simplified system may be used to make predictions on the design of a tissue-engineering scaffold and the appropriate operating regime that ensures a desired level of tissue growth.


Introduction and motivation
In vitro tissue engineering aims to create functional tissue and organ samples external to the body to replace damaged or diseased tissues and organs needed in multiple clinical therapies (Rose and Oreffo 2002;Martin et al. 2004). By using autolo-gous cell sources, often seeded within a porous scaffold that acts as a template for the developing tissue, tissue-engineered products have many advantages over traditional approaches such as donor tissue and organ transplantation that can elicit an adverse immune response. The development of the growing tissue construct-the combination of scaffold, cells and extracellular matrix (ECM) with biological cues-often occurs within a bioreactor. In perfusion-bioreactor systems, a flow of culture medium (nutrient-rich fluid) is driven through the porous scaffold. This serves two purposes. Firstly, the flow imposes fluid mechanical forces (shear stress, pressure) on mechanosensitive cells found in, for example, bone, cartilage, muscle, liver and blood vessels (Mullender et al. 2004). The mechanical environment that cells experience affects their differentiation, proliferation, orientation, gene expression and a host of other activities: the mechanical stimuli are integrated into the cellular response via a process known as mechanotransduction. For example, fluid shear stress has been shown to enhance extracellular matrix formation in vitro (Bakker et al. 2004;Klein-Nulend et al. 1995;You et al. 2000), and stimulation by hydrostatic pressure promotes stem cell differentiation down the chondrogenic lineage to form cartilaginous tissues (Elder and Athanasiou 2009). Additionally, fluid flows are employed to ensure adequate solute delivery to the cells within the porous scaffold. A variety of dynamic flow-perfusion bioreactors have been developed, including direct and free-flow perfusion systems, and the rotary cell culture system (RCCS) (see Martin et al. (2004) and references therein). By selecting an appropriate bioreactor system, tissue engineers aim to provide the optimum stimulatory environment for in vitro tissue growth to engineer a tissue construct that remains functional for significant periods of time.
In addition to the selection of perfusion bioreactor, the scaffold properties are of central importance to the tissue engineering approach. Scaffolds with a highly porous, interconnected structure encourage cell penetration, vascularization of the construct from the surrounding tissue (when implanted in vivo) and efficient mass transfer of nutrients and waste products. The scaffold material must be compatible with the host material so that it does not elicit an adverse immune response upon implantation (Salgado et al. 2004). The scaffold plays the role of the extracellular matrix (ECM) in tissues and defines its mechanical properties (e.g. the collagen and elastin fibres present in ECM, or the calcium hydroxyapatite that lends bone tissue its rigidity). One challenge is to select the appropriate material properties for the scaffold (porous, poroelastic, viscoelastic, etc.) to ensure that the cells experience the appropriate mechanical environment when the construct is subjected to an applied load. Another challenge is to ensure that the scaffold is sufficiently porous without compromising its mechanical integrity. Additionally, the rates of scaffold degradation (e.g. due to hydrolysis) must be matched to the rate of nascent tissue growth in vivo so that once the construct is in vivo, the artificial scaffold is replaced by native ECM and the mechanical integrity of the construct is maintained. Finally, both topographical and biochemical characteristics of the porous scaffold have been exploited to enhance cell-scaffold adherence, and to direct cell movement or differentiation (Weiss 1945;Salgado et al. 2004): for example, scaffolds can deliver growth factors or DNA (Sipe 2002), or contain specific cellular recognition molecules (Freed and Vunjak-Novakovic 1998). For a comprehensive discussion of scaffolds for tissue engineering applications see Atala et al. (1997), Hutmacher (2006), Hollister (2005 and references therein. The motivation for this work is the experimental observation that tissue growth is strongly affected by the geometrical features of the pores in an artificial 3D tissue scaffold (Rumpler et al. 2008), which has formed the focus of many recent studies (see, for example, Bidan et al. 2013;Werner et al. 2017;Lo et al. 2016;Kommareddy et al. 2010, among many others). For example, Bidan et al. (2013) emphasized how geometric constraints may guide tissue formation in vitro and showed that optimizing scaffold architectures may improve tissue formation independently of the scaffold material used. Furthermore, Werner et al. (2017) systematically investigated the influence of 3D surface curvature on cell and nucleus morphology and the subsequent effect on the migratory and differentiation behaviour of human mesenchymal stem cells. Key to successful tissue engineering is the coordination of the various cellular cues to enhance tissue growth. Mechanistic mathematical modelling has an important role to play in this regard, elucidating mechanisms and providing quantitative information about the cellular microenvironment (e.g. the shear-stress field) that is not straightforward (or even possible) to measure experimentally. Mathematical models may be validated against measurable experimental data, such as perfusion flow rate or bioreactor outlet nutrient concentration, and then exploited to reveal details of the mechanical and biochemical fields within the scaffold and bioreactor system. The quantitative models can be used to predict the outcome of a particular experimental scenario (limiting the need for numerous and expensive bioreactor experiments, potentially saving time and resources) and to optimize tissueengineering protocols.
There has been a large body of work using mathematical models to understand the relationship between flow, associated transport and tissue growth in perfusion-bioreactor systems with a porous scaffold. Numerous macroscale continuum models for bioactive porous media have been proposed for artificial tissue growth, with many accounting for the influence of fluid shear stress on tissue growth (see O'Dea et al. (2012) and references therein). 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 appropriate constitutive laws are specified to capture the interactions between the phases, enabling a wide variety of biological systems to be modelled (O'Dea et al. 2010;Pearson et al. 2014). An alternative approach is to include the effect of cell growth on the fluid flow implicitly, by prescribing the scaffold porosity to be a function of cell density, which satisfies a conservation-of-mass equation (Coletti et al. 2006;Shakeel et al. 2013). Fluid flow through the porous scaffold is then governed by Darcy's law. In these macroscale approaches, constitutive assumptions relate microscale processes, such as cell growth, to macroscale parameters. Alternatively, homogenization techniques have been proposed to capture the precise relationship between macroscopic parameters in continuum models, such as porosity, and porescale growth and flow processes (Shipley et al. 2009;Penta et al. 2014;O'Dea et al. 2015). As an alternative to continuum models, network models have been proposed for artificial tissue growth, particularly where the scaffold biomaterial contains relatively few pores (Krause et al. 2017a, b). These continuum and network approaches complement com-putational and algorithmic approaches employed to model various aspects of tissue engineering; see Geris (2013) for a comprehensive review.
The above studies consider tissue growth on the scale of the tissue construct. Although there are studies that have mapped the shear stress at the pore scale for realistic flows through a perfusion bioreactor (see, e.g. Raimondi et al. 2005), to the best of our knowledge, no mathematical studies have considered the interplay between fluid shear stresses and pore curvature on cell growth at the scale of an individual pore. In this paper, we will focus on the combined roles of fluid shear stress and topographical cues, in particular pore curvature, on tissue growth within a single pore of a scaffold. The paper is laid out as follows: in Sect. 2, we introduce a mathematical model for flow through a pore of a tissueengineering scaffold. The governing equations describing the flow of nutrient-rich solution and tissue growth are presented in Sects. 2.1 and 2.2, respectively. The model is analyzed in Sect. 3 by using an asymptotic approach exploiting the small aspect ratio of the pore. Asymptotic and numerical techniques are used to obtain solutions for the tissue growth in Sect. 4. We also discuss which initial internal pore morphology gives more tissue growth in Sect. 4.4. Finally, we conclude in Sect. 5 with a discussion of our model and results in the context of pores within a real tissue-engineering scaffold.

A mathematical model for flow within a scaffold pore
We consider a single scaffold pore lined with cells, through which nutrient-rich culture medium flows. Over time, the cells proliferate and the pore will constrict as the tissue layer lining it thickens. We consider the dynamics of the system on the timescale of cell proliferation, as the cells lining the pore wall grow on a much longer timescale than that associated with the transport of fluid through the pore (see values presented in Table 1 later). Thus, the pore wall may be considered stationary when considering the fluid flow, a quasi-static assumption.
The fluid-cell-layer interface is of lengthL and has typical radiusR. (Hats will be used throughout to distinguish dimensional quantities from their dimensionless equivalents.) We define the aspect ratio =R/L 1. We will explore how modifications in fluid-cell-layer interface geometry influence the resulting tissue growth. We will assume that the initial cell layer is sufficiently thin that its internal configuration reflects the underlying substrate geometry, which we assume to be approximately circular with small azimuthal and axial variation. We consider the regime in which the fluid-cell-layer interface remains nearly circularly cylindri-cal so that deviations of the circumference geometry from circular are small (see Fig. 1).
Nutrient is assumed to be present in excess so that no significant depletion occurs over the length of the pore. In many applications, the flux of fluid is held constant to ensure a continued supply of nutrient (Pearson et al. 2014;Shakeel et al. 2013); we assume this throughout. Specifically, we prescribe the inlet fluxQ i and the downstream pressureP d . The nutrient-rich culture medium is modelled as an incompressible Newtonian viscous fluid of (dynamic) viscosityμ and densityρ.
We work in cylindrical coordinates (r , θ,ẑ), whereẑ is aligned with the pore axis, and e i correspond to the unit vectors in the i direction. The fluid-cell-layer interface is described byr =â(θ,ẑ,t), and the initial configurationr = a(θ,ẑ, 0) is prescribed.

Fluid dynamics
The fluid flow has velocityû =ûer +ve θ +ŵeˆz and pressurê p. Neglecting flow inertia within the pore (Chung et al. 2007;Lemon et al. 2006), we assume the flow is governed by the Stokes equations. We introduce the following scalings: The dimensionless governing flow equations are then holding in the flow domain 0 ≤ r ≤ a(θ, z, t), 0 ≤ z ≤ 1. These equations must be solved subject to the following boundary conditions: enforcing no slip and no penetration at the fluid-cell-layer interface, and the symmetry conditions We assume that the (dimensionless) pressure drop across the length of the pore is given by ζ(t), which will increase monotonically to sustain the constant flux as the pore constricts under cell growth. The system is thus closed by applying the boundary conditions and enforcing constant fluid influx, 2π 0 a 0 w(r , θ, 0, t)r dr dθ = π.
Note that setting ζ(t) = 1 in (8), and dropping the condition (9), describes the alternative scenario in which the pressure drop is constant.

Tissue growth
For the cell proliferation (tissue growth) within the pore, we prescribe a phenomenological law based on experimental observations. Specifically, cells proliferate more quickly when exposed to higher shear stresses at their surface (mechanotransduction) and in regions where the underlying substrate has higher curvature (see, for example, O'Dea et al. (2010) and references therein; and Rumpler et al. (2008)). In the original dimensional variables, we propose Here,κ =∇·n is the mean curvature, with∇ the dimensional gradient operator and n the unit normal to the fluid-cell-layer interface (in the direction pointing into the fluid), given by The characteristic growth rate (m 2 s −1 ) is given byλ, and the function f captures the specific dependence of tissue growth on the total fluid shear stress at the fluid-cell-layer interface, σŝ. We will explore physically realistic expressions for f in Sect. 4 and discuss the limitations of the growth law (10) in Sect. 5. Without the shear stress dependence ( f (σˆs) = constant), (10) is precisely the growth law considered by Rumpler et al. (2008). Note also that, for cylindrical pores, κ = 1/a and (10) then implies that the rate of change of tissue volume depends only on shear stress, i.e. Dimensionless fluid-cell-layer interface aspect ratio 0.05-0.5 All parameters depend on the specific application; here, sample values have been taken, andλ was extracted from the gradient of projected tissue area versus time graph in figure 4b of Rumpler et al. (2008) At the fluid-cell-layer surface, the stress vectorτ is given bŷ wherê is the Cauchy stress tensor for the fluid, where I is the identity matrix. The magnitude of the stress vector in the normal direction is given bŷ The total shear stress at the fluid-cell-layer interface,σˆs, is then given bŷ Using the non-dimensionalization presented in (1) together with the additional scaleŝ gives the dimensionless growth law The unit normal to the fluid-cell-layer interface, using (11), is given in dimensionless form by 3 Model analysis

Asymptotic ansätze
The full system outlined in Sect. 2 can be solved numerically for a specified growth function f . However, this would be computationally costly, limiting the parameter ranges that can be explored and affording only limited insight into the effects of model parameters. We therefore take an alternative asymptotic approach, exploiting the smallness of the typical pore aspect ratio, =R/L 1. This enables us to make analytical progress, leading to a more tractable model, which can be solved cheaply to provide insights into the dynamics of the system.
In addition to the assumption of small pore aspect ratio, we assume that the initial cell layer is sufficiently thin that its internal configuration reflects the underlying substrate geometry, which we assume to be approximately circular with small azimuthal and axial variation. This motivates relegation of spatial variations in the fluid-cell-layer interface to first (axial) and second (azimuthal and axial) order in the aspect ratio. Below we find that the leading-order fluid velocity is unidirectional and independent of the axial coordinate z. As a result, the leading-order shear stress at the fluid-cell-layer interface is independent of z. We make the quasi-static assumption of slow tissue growth compared with the timescale of nutrient transport through the scaffold, and the further assumption that there is a plentiful supply of nutrient so that all parts of the growing tissue experience the same concentration. As a result, the cell layer grows in a spatially uniform manner with time to leading order in the pore aspect ratio. Continuing the analysis to higher order in the aspect ratio, we determine the dependence of the tissue growth on the axial and azimuthal coordinates at first and second order in the pore aspect ratio, respectively.
We therefore assume that the fluid-cell-layer interface may be expressed as We further assume that a 2 is of the form where Λ 2 (z, t) and Υ 2 (z, t) are functions to be determined and n is an integer. This form corresponds to periodic perturbations of the pore from a circular cross section that has n lobes. While this ansatz clearly imposes constraints on the type of tissue growth we study, we will find that this strikes a balance between allowing us to make substantial analytical progress while retaining a suitable level of generality to the form of the cross-sectional profile. We assume that the initial fluid-cell-layer interface shape a(θ, z, 0) is given, then determine its subsequent evolution via Eq. (18).
In the following section, we pose asymptotic expansions of the form and similarly for v, w, p, ζ , σ s and κ, and systematically examine the system at each order to determine reduced equations.

Leading-order analysis
For the fluid-cell-layer interface configuration (20), the flow will be independent of θ to leading order in . Substituting the asymptotic expansions (22) into Eqs.
(2)-(5) and conditions (6)-(9) and retaining leading-order terms give a system of equations for the leading-order velocities, u 0 , v 0 and w 0 , pressure p 0 , and ζ 0 (the leading-order pressure at the pore inlet, determined as part of the solution in our specified-flux scenario):

First-order analysis
At O( ), based on the ansatz (20), we seek a θ -independent solution to the governing equations, finding

Second-order analysis
Considering Eq. (4) at O( 2 ) motivates seeking a solution w 2 of the form Using (4) and (6), we find Imposing the boundary conditions (6)-(8) and the flux condition (9) gives, after some manipulation, and the axial velocity is given by Note that u 2 and v 2 may also be determined, but are not required in the analysis that follows.

Asymptotic analysis for the tissue growth
Equation (19) yields the asymptotic expansion of the unit normal, Substituting (17), (32) and (33) into (16), we obtain the dimensionless shear stress at the fluid-cell-layer interface as where The dimensionless mean curvature is given by where We integrate the growth-law Eq. (18), to yield an explicit expression for the fluid-cell-layer interface shape, Substituting our asymptotic expansions for a, σ s and κ, Eqs. (20), (34) and (36), respectively, into (38), give Here, primes denote differentiation with respect to the argument and arise due to a Taylor expansion of f , while σ s0 , σ s1 , σ s2a , σ s2b , κ 0 and κ 2 are given by (35) and (37).

Results
In this section, we present simulations of the model (39)-(42) for a physically realistic growth law. Specifically, we choose f (see Eq. (18)) to be of the form: where F 1 < F 2 , σ 1 < σ 2 and m 1. Based on the shearstress scaling in (17), σ 1 and σ 2 will typically be order-one quantities (Rumpler et al. 2008). For large values of m, this corresponds to a function that has three regions in which it is approximately constant, connected by rapid but smooth transition regions (see Fig. 2). Similar constitutive laws have been considered elsewhere (see, for example, Krause et al. 2017a, b;Pearson et al. 2014). This captures low growth at small shear stresses ( f ≈ F 1 for σ < σ 1 ), enhanced growth at intermediate shear stress ( f ≈ F 2 > F 1 for σ 1 < σ < σ 2 ), and the adverse impact of high shear stress on growth ( f ≈ 0 for σ > σ 2 ). We note that our analysis readily generalizes for nonzero growth, or even negative growth (capturing cell death and degradation), for high stresses, but we do not consider this here. We note that as m becomes large, gradients in f become large (order m). However, our Taylor series expansion in (39)-(42) remains valid, since the window over which the gradient is large is only of width 1/m. Thus, our asymptotic expansion in holds, irrespective of the size of m.
In Case (i), tissue growth begins in the first (low growthrate) phase at a rate F 1 . The tissue growth leads to pore shrinkage, hence (with the imposed constant fluid flux), increased shear stress. The growth rate thus transitions into the second phase, where growth occurs at a higher rate F 2 , before finally transitioning into the third phase, where growth ceases. In Case (ii), tissue growth starts in the second phase (at a rate F 2 ) and, as shear stress increases, moves into the third phase, where growth stops. Note that for most of our simulations, we choose the values F 1 = 1, F 2 = 3, σ 1 = 7 and σ 2 = 15 for the growth function f given by (43). How the choice of these parameters affects model results is discussed further below. Our model contains two dimensionless geometrical parameters: the pore aspect ratio, , and the azimuthal wavenumber n that appears in the perturbation a 2 in Eq. (21). The pore aspect ratio, , is unknown for the experiments of Rumpler et al. (2008) and varies between scaffolds: here, we set = 0.2 for most of our simulations to highlight more clearly the effects of pore shape on the tissue growth in the θ direction. This appears to be a reasonable choice, based on the qualitative comparison made in Fig. 7 later. We consider a range of positive integer values for n to gain insight into the role it plays in the total amount of tissue growth and final patterns obtained.

Case (i): s0 (0) < 1 < 2
While the system (39)-(42) may be readily solved numerically for any suitably smooth differentiable function f , we begin by considering the specific functional form of f , given by (43), when m → ∞. This corresponds to the case where the tissue growth rates undergo jump discontinuities at critical values of the shear stress. We can then solve the system analytically for the three separate constant-growthrate regions corresponding to σ s < σ 1 (when f = F 1 ), σ 1 < σ s < σ 2 (when f = F 2 ), and σ s > σ 2 (when f = 0). Finally, we connect these regions to determine the full solu-tion. This analytical approach provides physical insight into the parametric dependence of the system growth.
From equations (23) and (35), σ s0 = 4/a 3 0 ; therefore, (37), (39) and (43) give when σ s0 (0) < σ 1 . Here, t * 01 is the time at which σ s0 = σ 1 , determined by a 0 (t * 01 ) = (4/σ 1 ) 1/3 , and so may be explicitly written as At this time, the cell growth transitions into the second phase, where σ 1 ≤ σ s0 < σ 2 . Equation (39) then gives where a 0 (t * 02 ) = (4/σ 2 ) 1/3 and Equation (47) provides an explicit prediction for the time at which tissue growth stops and its dependence on the model parameters σ 1 , σ 2 , F 1 and F 2 in the limit m → ∞. We note that, when we substitute for t * 01 using (45), there is a rather complex dependence on each of the parameters that appears in our tissue-growth ansatz, (43). However, importantly, order-one changes in any of these parameters will lead to order-one changes in the time taken for the tissue growth. This highlights the need to determine the value of these parameters for a given experimental scenario if this model is to form a predictive tool. We emphasize that in the subsequent examples, the values taken for these parameters are illustrative, but give qualitatively similar behaviour to results in the literature. Repeating the previous procedure for t ≥ t * 02 , we obtain the leading-order fluid-cell-layer interface position: We thus see that, to leading order, the position of the fluidcell-layer interface is (4/σ 2 ) 1/3 . This means that, for a given underlying pore structure, the total amount of tissue grown is independent of F 1 , F 2 and σ 1 . These additional parameters only control the rate at which the final state is approached.
Note that when substituting the functional form of f defined in (43) in Eqs. (40) and (42) we find that a 1 and Υ 2 do not depend on time and therefore retain their initial values throughout, i.e. a 1 (z, t) = a 1 (z, 0) and Λ 2 (z, t) = Λ 2 (z, 0). (This is because the integrals in (40) and (42) involve first and higher derivatives of the growth function f , which are thus zero.) This allows for the possibility of axial variations in the equilibrium fluid-cell-layer interface. Following a similar process to that used above to determine a 0 , Λ 2 (z, t) may be found from (35), (37), (41) and (43) as Thus, the evolving cross-sectional profile of the tissue structure is given by where a 0 (t) and Λ 2 (z, t) are given by (48) and (49), respectively. Note that if the nutrient solution is supplied under conditions of specified pressure drop rather than specified flux (so that ζ(t) = 1), then an identical analysis holds, except with We now compare the above (rather general) analytical solution with numerical simulations to the system (39)-(42) in the appropriate parameter regimes in the following analysis. We relax our assumption of m → ∞ in the expression for the growth function f given in (43) and consider finite values of m so that the tissue growth transitions smoothly between the three growth rate regions as shear stress increases. We first consider an initial fluid-cell-layer interface that is cylindrical to leading order, but for which higher-order corrections are all chosen to be linearly decreasing along the pore axis (see (31)): We choose σ 1 such that the initial leading-order shear stress σ s0 (0) < σ 1 , and so, as the fluid-cell-layer grows, with atten-dant increase in shear stress, we sweep through all three regions of the growth function f , given by (43).
In the simulation illustrated in Fig. 3, the wavenumber n = 4 is used in the azimuthal perturbation to the fluid-cell-layer interfacial profile (see (21)) and the initial fluid-cell-layer interface configuration is given by Eqs. (52). The black curves in Fig. 3a show the numerical solution to (39) for the leading-order fluid-cell-layer interface a 0 (t) versus time t, for several different values of m, while the red graph shows the analytical solution (see (48)). Our results here confirm that as the value of m increases, the numerical solution to (39) converges to the analytical solution for a 0 in (48). Furthermore, Fig. 3b confirms that when m 1 (m = 1000 is used here), the numerical solution to (41) for Λ 2 is almost identical to the analytical solution in (49). Therefore, for the remainder of the paper, we use the analytical solutions (m → ∞) to (39)-(42) for both Case (i) and (ii). Figure 4 illustrates the analytical results to (39)-(42) (m → ∞) as obtained in (44)-(50) for the initial fluid-celllayer interface configuration in (52). (The results of Fig. 4a, c were already shown in Fig. 3, but we now discuss them in more detail.) The leading-order radius of the fluid-celllayer interface, a 0 (t), decreases approximately linearly with t (Fig. 4a), which is consistent with a Taylor expansion of (44) for small time. The leading-order shear stress, σ s0 , also increases approximately linearly, until it surpasses the value σ 1 and we enter the regime of maximal tissue growth rate (σ 1 ≤ σ s0 < σ 2 ) (Fig. 4e). The fluid-cell-layer interfacial radius then decreases rapidly so that the shear stress also rapidly increases further, until finally it exceeds σ 2 and the growth is halted (see the solution for a 0 (t) in Fig. 4a).
Since the function f , given by (43), is piecewise constant for our chosen limit m → ∞, a 1 and Υ 2 do not vary in time, as discussed in our analytical work above (see Fig. 4b, d). The evolution of Λ 2 (z, t) is shown in Fig. 4c (see (49)): this function retains its linear dependence on the axial coordinate z, while approaching zero smoothly as time progresses. The evolution of the leading-and second-order shear-stress components, σ s0 and σ s1 , is shown in Fig. 4e, f, which are obtained from (23) and (35).
The overall picture may be visualized by considering the fluid-cell-layer interfacial profile at several different cross sections along the scaffold, both initially, and when the system has reached equilibrium (Fig. 4g, h, respectively). As observed, the tissue will grow so that ultimately the tube is circular in cross section, but the radius may vary along the axis at O( ) due to the initial choices of a 1 and Υ 2 .
We next consider an initial configuration in which the initial profiles for a 1 , Υ and Λ all increase rather than decrease with distance along the pore axis (the direction of flow):  (52) with a 0 (0) = 0.9, a 1 (z, 0) = −z −0.5, Λ 2 (z, 0) = −z +2, Υ 2 (z, 0) = −z +2 and growth function f given in (43) with F 1 = 1, F 2 = 3, σ 1 = 7, σ 2 = 15. The black graphs in a show the numerical solution for a 0 for several dif-ferent values of m, while the analytical solution (m → ∞) is given by the red dashed graphs (see (48)). b Λ 2 (z, t) at several different times, where black and red graphs are numerical (for m = 1000) and analytical (m → ∞) solutions, respectively. Here, t f = 0.25 is chosen to be sufficiently large that the system has reached an approximate steady state and tissue growth stops a 0 (0) = 0.9, (53a) (53d) As in the previous example, during the tissue growth, the shear stress attains values in all three regions of the domain of the growth function f , taking the tissue growth through regimes of moderate, fast, and ultimately cessation of growth.
The final fluid-cell-layer interfacial profile is again circular in cross section, though with radius that varies along the length of the pore. The evolution towards this final configuration is shown in Fig. 5.

Case (ii) 1 < s0 (0) < 2
When σ 1 < σ s0 (0) < σ 2 , cell growth begins in the second of the three regions of the domain of the growth function f illustrated in Fig. 2 (fast initial growth). As in Sect. 4.1, we examine the case m → ∞ in which analytical progress is possible. Substituting (23), (35), (37) and (43) into (39)-(42) and following a similar approach to that of Sect. 4.1 gives and where Note that a 1 and Υ 2 again do not depend on time, i.e. a 1 (z, t) = a 1 (z, 0) and Λ 2 (z, t) = Λ 2 (z, 0), as discussed earlier in Sect. 4.1. As in Case (i), the position of the fluidcell-layer interface is given, to leading order, by (4/σ 2 ) 1/3 , (54). Hence, for a given underlying pore structure, the leading-order amount of tissue grown will be independent of F 1 , F 2 and σ 1 . Equation (56) provides a prediction for the time at which tissue growth stops in this regime and its dependence on the model parameters in the limit m → ∞. Since this scenario is simpler than Case (i), the dependence of the time taken for the tissue to complete growth is less complex. In particular, in this case the time to cessation of growth is independent of σ 1 and F 1 , and depends on F 2 in an inversely proportional fashion. Nevertheless, as outlined in Case (i) in Sect. 4.1, there is a clear need to determine the values for σ 2 and F 2 if this model is to form a reliable predictive tool.  (44)-(50) for Case (i), Sect. 4.1 σ s0 (0) < σ 1 < σ 2 : a-d are profiles for a 0 , a 1 , Λ 2 and Υ 2 , respectively, for the fluid-cell-layer interface (20), with = 0.2, n = 4 and initial conditions (52), with a 0 (0) = 0.9, a 1 (z, 0) = −z − 0.5, Λ 2 (z, 0) = −z + 2 and Υ 2 (z, 0) = −z + 2; e, f are shear stresses at leading and orders, σ s0 (t) and σ s1 (t) (see (23) and (35)), respectively, and g, h are the fluid-cell-layer interfacial cross sections at initial and final times (t f = 0.25, chosen to be sufficiently large that the system has reached an approximate steady state and tissue growth stops), respectively, for the growth function f given in (43) with F 1 = 1, F 2 = 3, σ 1 = 7, σ 2 = 15 and m → ∞. Here, (x, y) = (r sin θ, r cos θ). The red dashed lines in a, e correspond to the transition times and the values of σ 1 , σ 2 , respectively In Fig. 6, we consider an initial configuration identical to that used in Fig. 4, given by Eq. (52), except that the leadingorder fluid-cell-layer interfacial radius is slightly smaller, a 0 = 0.8. This ensures that the initial leading-order shear stress σ s0 (0) is sufficiently large to lie within the second region of the domain of the growth function f in (43) for identical values of σ 1 , σ 2 , F 1 and F 2 . The tissue growth here is initially faster than in the previous cases: the fluid-celllayer interfacial profile shrinks rapidly and we quickly transit  (44)-(50) for Case (i), Sect. 4.1 σ s0 (0) < σ 1 < σ 2 : a-d are profiles for a 0 , a 1 , Λ 2 and Υ 2 , respectively, in the fluid-cell-layer interface (20), with = 0.2, n = 4 and initial conditions (53) with a 0 (0) = 0.9, a 1 (z, 0) = z − 0.5, Λ 2 (z, 0) = z − 2 and Υ 2 (z, 0) = z − 2; e, f are shear stresses at leading and orders, σ s0 (t) and σ s1 (t) (see (23) and (35)) respectively, and g, h are fluid-cell-layer interface cross sections at initial and final times (t f = 0.25, chosen to be sufficiently large that the system has reached an approximate steady state and tissue growth stops), respectively, for the growth function f given in (43) with F 1 = 1, F 2 = 3, σ 1 = 7, σ 2 = 15 and m → ∞. Here, (x, y) = (r sin θ, r cos θ). The red dashed lines in figures a, e correspond to the transition times and the values of σ 1 , σ 2 , respectively to region three of the domain of the growth function, in which growth is suppressed.
Note that, similar to what was shown in Fig. 3, excellent agreement was again obtained between the numerical solutions of the system (39)-(42) (for m = 1000) and the analytical asymptotic prediction (54) and (55) for the leading-order fluid-cell-layer interfacial radius a 0 as a function of time and the perturbation function Λ 2 (z, t) versus z at several different times (results not shown here).  (44)-(50) for Case (ii), Sect. 4.2, σ 1 < σ s0 < σ 2 : a-d are profiles for a 0 , a 1 , Λ 2 and Υ 2 , respectively, in the fluid-cell-layer interface (20), with = 0.2, n = 4 and initial conditions (52) with a 0 (0) = 0.8, a 1 (z, 0) = −z − 0.5, Λ 2 (z, 0) = −z + 2 and Υ 2 (z, 0) = −z + 2; e, f are shear stresses at leading and orders, σ s0 (t) and σ s1 (t) (see (23) and (35)), respectively, and g, h are fluid-cell-layer interface cross sections at initial and final times (t f = 0.25, chosen to be sufficiently large that the system has reached an approximate steady state and tissue growth stops), respectively, for the growth function f given in (43) with F 1 = 1, F 2 = 3, σ 1 = 7,  Figure 7 shows a direct comparison between experimental results visualized under a confocal laser scanning micro-scope by Rumpler et al. (2008), and analytical solutions of our model for several different values of the polygonal cross section n (the wavenumber in the azimuthal perturbation to the initial fluid-cell-layer interfacial profile (see (21)). The results show that our model gives a reasonable qualitative agreement with the experimental results, illustrating the tendency of the fluid-cell-layer interface to become ultimately circular, independent of the underlying pore geometry. Note that here the experimental tissue grown (in Fig. 7a-c) is shown after 21 days, while in our simulation the dimensionless time at which the tissue growth is observed to stop is t f = 0.25 (see , which corresponds to a dimensional time of 20.5 days using the time scale identified in Table 1. Furthermore, as shown in Fig. 7 (and as expected from the form of the growth law), more tissue grows at the corners of the scaffold structure so that the fluid-celllayer interface subsequently becomes more circular. This behaviour leads to the formation of a final round fluid-celllayer interface, regardless of the original shape (at least for the parameters chosen here).

Optimization and parameter sensitivity
A question of interest to tissue engineers is the optimal scaffold configuration, in terms of internal pore morphology (pore size and shape) that maximizes net tissue yield over the course of an experiment. To answer this question, we first identify a metric for tissue growth. The total volume V (t) of tissue grown within the pore by time t is given by (57) To gain insight into the influence of the pore shape on the tissue growth, we first consider the dependence of V (t) on the value of n, the number of azimuthal modes in the internal fluid-cell-layer configuration within the pore, in Fig. 8a.
These results indicate that introducing regions of high curvature in the fluid-cell-layer interfacial structure enables tissue to be formed more quickly in the structure over the intermediate growth period. During the initial growth phase, V (t) is an increasing function of n: for t = 0.025, for example, the volume of tissue with n = 6 is approximately double that with n = 2. However, as t increases, the volumes for all n values equalize, and the final amount of tissue growth is largely independent of the number of lobes, n. Thus, we may conclude that if the rate of tissue growth is a primary objective, then increasing the number of lobes in the tissue-engineering construct is advantageous. However, if the concern is only with the final volume of tissue grown, then the pore shape is less important.
In order to further study the model's parameter sensitivity, we also investigate how our results depend on , the fluidcell-layer interfacial aspect ratio. A graph showing the total amount of tissue grown within the pore, V (t), as a function of , is shown in Fig. 8b. Here, our results show that the final amount of tissue grown increases with (though only by about 10% over the model's range of validity). The variation of total tissue grown with manifests itself through Λ 2 , given by (49) and (55) for cases (i) and (ii), respectively.
The temporal dependence of tissue growth on the parameters F 1 and F 2 , and σ 1 and σ 2 is shown in Fig. 8c, d, respectively. The final position of the fluid-cell-layer interface was seen to be given to leading order by (4/σ 2 ) 1/3 (Eqs. (48) and (54) for Cases (i) and (ii), respectively). This indicates that for a given underlying pore structure, the total amount of tissue grown is independent of F 1 , F 2 and σ 1 to leading order, which is confirmed in Fig. 8c, d.

Conclusions
We have presented a simplified mathematical model for the growth of tissue in a tissue-engineering scaffold to gain insight into the effect of pore morphology on the tissue growth. The flow of nutrient to the cells was captured by the Stokes equations, while the cell proliferation was modelled by a law that accounted for the effects of the underlying fluidcell-layer interface curvature, and fluid shear stress at the growing tissue surface, on the rate of tissue growth. Exploiting the geometrical features of a typical scaffold, namely a structure composed of a series of pores that are nearly cylindrical, allowed us to proceed via an asymptotic approach that led to a reduced system of four simple differential equations for the tissue growth. The resulting equations were analyzed numerically and analytically for a typical growth law, and an analytic expression was obtained when the growth law was piecewise constant. The analytic solution was shown to agree with the numerics, supporting its use for rapid testing and analysis of the system in comparisons with experimental data from Rumpler et al. (2008).
The number of azimuthal modes in the initial fluid-celllayer interface profile, n (which we take as a proxy for the shape of the underlying scaffold pore), was found to have a minimal effect on the final amount of tissue grown. These results suggest that focusing effort into engineering the pore geometry may be beneficial if the rate of tissue growth is a principal concern, but if maximizing the final volume of tissue that can be grown is the main objective, then the details of the shape of the (nearly cylindrical) pore within the scaffold are not important. In contrast, the total amount of tissue grown depends on , the fluid-cell-layer interface aspect ratio, with the amount of final tissue grown increasing with increasing .
The results of this work, once calibrated against available data, offer a simple framework for testing the behaviour of different scaffold pore geometries without the need for many complex experiments to be conducted. The model should ultimately allow for the prediction of regimes in which tissue growth may be enhanced, offering improvements and  Rumpler et al. (2008). d-f Upstream side (z = 0) of the fluid-cell-layer interfacial cross section at initial and final times (t f = 0.25, chosen to be suitably large such that the system has reached an approximate steady state) for n = 3, n = 4 and n = 6, respectively, with = 0.2 and for the initial conditions in (52) with a 0 (0) = 0.9, a 1 (z, 0) = z − 0.5, Λ 2 (z, 0) = z − 2 and Υ 2 (z, 0) = z − 2, the growth function f given in (43) with F 1 = 1, F 2 = 3, σ 1 = 7, σ 2 = 15 and m → ∞. Here, (x, y) = (r sin θ, r cos θ) additional insight into the current operating strategies used. To compare our theoretical results with experiments requires an estimate ofλ, the growth-rate parameter in Eq. (10), for the dataset under consideration. This parameter depends on many experimental details, including cell type, the physical characteristics of the tissue-engineering scaffold and the concentration of the nutrient in the culture fluid. In principle, λ could be measured or estimated directly from a suitable experimental dataset, as we did for the data of Rumpler et al. (2008). Here,λ was estimated by extracting the gradient of the projected tissue area versus time graph in Fig. 4(b) of Rumpler et al. (2008). The time taken for the tissue to form and the total amount of tissue grown were inherently linked to the parameters in our tissue-growth ansatz. For this model to form a reliable predictive tool, the parameters that appear in the growth ansatz need to be determined experimentally. While this is possible in principle by examining temporal data, in this paper we had access only to final-time snapshots. A close collaboration between experimentalists would be necessary to pin down these parameters more accurately, and represents a future research direction.
The model presented in this paper has the advantage of simplicity, but obviously fails to capture many details of the real experimental system, which may be important in some scenarios. For example, throughout this study, the nutrient is assumed to be present in excess so that cells always have sufficient nutrient to proliferate, and no significant nutrient depletion occurs over the length of the pore. If nutrient is present in low concentration, this may not be true. In such cases, our model could be generalized to capture the effect of nutrient concentration gradients Cummings 2017, 2018). While significant gradients in nutrient concentration are unlikely to arise at the single-pore scale, they could certainly be present over the scale of an entire scaffold, consisting of many pores. One could consider scaling-up our model to describe an entire scaffold, with the nutrient flux varying in both time and spatial location within the scaffold. The growth parameterλ [see (10)] would then be considered a function of the local nutrient concentration, again informed by available data.
Furthermore, our growth model is empirical, based purely on experimental observations of how a tissue interface grows   (57)) versus time, for fluid-cell-layer interfacial profile (20), with the initial fluid-cellinterface configuration given in (52), a 0 (0) = 0.9, a 1 (z, 0) = −z −0.5, Λ 2 (z, 0) = −z + 2, Υ 2 (z, 0) = −z + 2 and for a several different values of the polygonal cross section n (the wavenumber in the azimuthal perturbation to the fluid-cell-layer interfacial profile (see (21)); b fluidcell-layer interfacial aspect ratio ; c F 1 and F 2 (see (43)); and d shear stresses σ 1 and σ 2 (see (43)). In (a), = 0.5; in (b), n = 4 and the growth function f is as given in Fig. 2 with F 1 = 1, F 2 = 3, σ 1 = 7 and σ 2 = 15. In (c) and (d), = 0.5 and n = 4. In all cases m → ∞ within pores of differing curvatures. Though these observations were validated across a range of different pores (Rumpler et al. 2008), one could argue that a more complete model is needed, which accounts for the details at the cell scale and the underlying mechanisms that drive cells to proliferate. Much remains to be done in terms of building such a detailed predictive model, which could offer improved insight into how the structural details of an underlying scaffold impact tissue engineering outcomes. Nonetheless, we believe that this model, which is characterized through a simple yet descriptive growth law that may be connected to data through the growth parameterλ and can be solved with minimal computational effort, offers some utility in addressing scaffold design questions.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecomm ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.