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

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

Introduction

In vitro tissue engineering is a form of regenerative medicine which often involves seeding cells into a porous bio-engineered scaffold that allows nutrient transport, structural support and a means for cell signalling activity (Chan and Leong 2008). Subject to the correct environment and growth factors, the cells will develop into a functional construct that can be used to restore damaged tissues and improve concepts in pharmaceutical research such as experimental drug therapy (Jensen et al. 2018). Employing contributions from an assortment of scientific fields, tissue engineering is considered an interdisciplinary practice that has the potential to benefit a substantial proportion of the global population with devastating soft tissue, bone and whole organ diseases (Dzobo et al. 2018). The field of tissue engineering has enjoyed many successes; for example, the generation, replacement and longevity of engineered bones and bronchial tubes derived from the recipients’ cells (Sato et al. 2008; Petite et al. 2000; Schimming and Schmelzeisen 2010). However, a shortage in the supply of donor tissue creates a demand on the field to make engineered tissue routinely clinically available (Levitt 2015). Whilst the field is rich in both theoretical and experimental knowledge, a lack of understanding regarding the processes by which cells assemble into tissues means that viable replacement constructs are only available in a minority of cases.

Extensive mathematical research has been undertaken to make sense of the complicated mechanisms within tissue growth. Some authors adopt a microscale approach, which can take the form of cellular automaton systems (Vitvitsky 2014; Lehotzky and Zupanc 2019; Youssef 2015) that seek to model interactions between a large number of individual cells. Whilst such systems can track the behaviour of cells, they can become computationally infeasible for tissue-scale simulations (Ermentrout and Edelstein-Keshet 1993). Some authors adopt a probabilistic approach (Fadai et al. 2019; Browning et al. 2019a, b) and exploit experimental data to explore ways in which tissue engineering techniques can be improved. For example, Sogutlu and Koc (2013) present a stochastic model to determine the expected required number of pores for each region of a tissue engineering scaffold.

Conversely, continuum models (O’Dea et al. 2010; Lemon et al. 2006; Breward et al. 2002; Eyles et al. 2019; Byrne et al. 2003 ) track the evolution of tissue constituents by employing systems of partial differential equations. Whilst continuum models cannot impose rules on individual cells, they can be derived by imposing averaging techniques (Drew 1983) on equations that govern cellular behaviour at a microscopic level, from which relevant mathematical techniques may be exploited to determine relationships between parameters and model behaviours. The study of in vitro tissue growth via continuum models is extensive, see O’Dea et al. (2012) and Klika et al. (2016) for reviews; of particular relevance to this study, however, Lemon et al. (2006) considers a continuum multiphase model to investigate how mechanical pressures within growing tissue influence the aggregation or dispersion of cells in a scaffold, and relates the existence of these regimes to the governing parameters. Lemon and King (2007a) examines travelling-wave solutions of the multiphase model formulated in Lemon et al. (2006) and find that in certain limits, the tissue propagates through the scaffold at a constant speed as either a forward or backward travelling wave, dependent on parameter values.

In this paper, we develop and analyse a continuum multiphase model that represents the development of tissue in vitro in an artificial scaffold. In our model, we aim to capture key features of tissue growth and extinction whilst developing a tractable formulation. In particular, we consider a porous flow description comprising a tissue cell phase, extra-cellular liquid phase and a scaffold phase, the former two being modelled as incompressible fluids and the latter as an inert solid. Tissue mechanics are accounted for by considering relevant constitutive assumptions in a similar fashion to those presented in Lemon et al. (2006) and Lemon and King (2007a). The model is reduced to a reaction–diffusion equation for the cell phase and a moving boundary condition for the tissue edge, after which travelling-wave, asymptotic and numerical methods are employed to deduce the resulting solution behaviour. The paper is constructed as follows. In Sect. 2, we formulate and subsequently reduce and non-dimensionalise the model. In Sect. 3, we present and discuss numerical solutions to the reduced model, which motivate the travelling-wave and asymptotic analyses conducted in Sects. 4 and 5. In Sect. 6, we draw some conclusions regarding the behaviour of the model and interpret the mathematical results in terms of the biological application.

Model Development

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

Governing Equations

We assume that cells proliferate and assemble daughter cells from the available liquid and that when cells die, they decompose and dissolve into the liquid phase. In view of these processes, it is reasonable to follow Lemon et al. (2006), Byrne et al. (2003), Breward et al. (2002) and Preziosi and Tosin (2003) (and many others) and assume the densities of n and w to be equal. Following these assumptions, the mass transfer equations can be represented aswhere is the net rate of cell proliferation. Adding the equations from (2) results in the overall conservation of mass conditionwhere (1) has been used to eliminate the time derivative and to replace w with

Noting that n and w are modelled as fluids and s as a porous scaffold, we take the interphase drags to be dominated by those with the scaffold and neglect that between the tissue and liquid. In view of this, we apply Darcy’s law to relate the velocity of the cell and liquid phases to their respective pressures. Following King and Franks (2004) and Eyles et al. (2019), we takewhere and represent the viscosity of the cell and liquid phases and K is the permeability of the scaffold.

Remaining consistent with Lemon et al. (2006) and Lemon and King (2007a, 2007b), we relate the cellular and extra-cellular liquid pressures viawhere represents extra pressures that arise due to cell–cell and cell–scaffold interactions. Since the scaffold is assumed to be inert and of uniform porosity, we suppress the dependence has on s from hereon for brevity. We note that combining (5) with the relations from (4) allows the elimination of and and provides

Initial and Boundary Conditions

Assuming the tissue to be symmetric about its centre (), we takeNaturally, the cell volume fraction is identically zero at the edge of the tissue:The moving boundary L(t) moves with the cell velocity, henceThe initial distribution of n and tissue boundary position, respectively, are denoted by

Model Reduction

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

Constitutive Assumptions

We now define constitutive assumptions for and that are suitable to describe tissue growth in a rigid scaffold. We assume that daughter cells are constructed via mitosis using the available liquid and that when cells die via apoptosis, they dissolve into the liquid. Thus, we takewhere and are the positive constant rates of cell mitosis and apoptosis, and (1) is used to replace w with Following Lemon et al. (2006) and Lemon and King (2007a, 2007b), an appropriate expression for isfor and positive constants and The first term in (16) represents repulsive forces exerted between the cells at high volume fractions, as characterised by the singularity at The second term represents the propensity for cells to disperse or aggregate, with taking a positive or negative value accordingly. The third term represents repulsive forces that occur due to cell–scaffold interactions, whilst the fourth describes attractive forces between the cells and scaffold. For simplicity, we take . We note that must be strictly positive to prevent negative diffusion in (12) and nonlinear degeneracy in (13). This is achieved when which is henceforth assumed. Physically, this corresponds to a tendency for cells to spread through the scaffold (Lemon et al. 2006).

Non-dimensionalisation and Parameter Values

We non-dimensionalise (12), (13) and the initial and boundary conditions from (8),  (10) and (14). By introducing the dimensionless variablesthe following dimensionless model results:where and . We also havewhere and In the proceeding, we dispense of the hat notation for clarity.

The parameter is shown in subsequent sections to be of crucial importance to the qualitative features of the model solutions. Physically, represents the difference between the scaffold porosity and the ratio between the cell death and growth rates. We note that the scaffold permeability parameter K,  as seen in (11) is not present in the non-dimensional model (18)–(20). However, the scalings selected in (17) imply that the dimensional tissue boundary position increases with the scaffold permeability. A linear stability analysis around the steady states of (18), and , is conducted in “Appendix A” and shows that the former is stable when and the latter when In view of this, we are primarily motivated to investigate (18)–(20) for different values of though variations in and will also be considered in part so as to deduce their optimal values for the stimulation of tissue growth. Unless otherwise stated, we take and we adopt the initial conditionsso that denotes the cell volume fraction at Following Lemon and King (2007a, 2007b), and unless otherwise stated, we set and the former corresponding to a scaffold with a porosity of 0.8 and is consistent with the experimental study presented in Malda et al. (2004).

Numerical Results

We present and discuss the numerical solutions for and L(t) from the PDE system (18)–(20), paying separate attention to the cases and For numerical convenience, we fix the moving boundary by introducing the variable transform so that which means (18) and (19) becomeSubject to the transformed boundary and initial conditions from (20), we numerically integrate (23) and (24) by discretising first- and second-order spatial derivatives using second-order finite differences. Upwind finite differences were used for the second term of (23). Temporal derivatives are numerically integrated by utilising ode23s in MATLAB.

For , as shown in Fig. 1a, b, we observe semi-infinite travelling waves in n and linear growth in L after a period of transient growth from their initial states. For , as shown in Fig. 1c, we observe n decaying from the initial data. Figure 1d shows unbounded growth in L. The inset shows L(t) and the function plotted against , from which we conclude that L grows logarithmically at large time. For , as shown in Fig. 1e, we observe that n decays from the initial data more quickly than for . The initial growth of L shown in Fig. 1f occurs due to the diffusion of n from the initial state; however, we observe the eventual formation of a steady state. Numerical simulations that are not included here suggest that travelling-wave and steady-state behaviour is exhibited by (18)–(20) for all and accordingly.

Fig. 1

Numerical solutions of the system of PDEs from (18)–(20), with a, c, e showing and b, d, f showing L(t) for different values of . The green lines in a, c, e represent the initial conditions from (22), whereas the black lines represent n for different values of t. The black arrows point in the direction of increasing time. The dashed red line in d is given by and highlights the logarithmic growth of L at large time. Parameter values: and (Color figure online)

Clearly, the case in which corresponds to effective tissue growth. This motivates a travelling-wave analysis of (18)–(20) for which is presented in Sect. 4, where we express the speed of the tissue edge in terms of the governing parameters. In Sect. 5, asymptotic solutions for n and L are found when so that the cell distribution and tissue speed are explicitly available. Whilst the case results in tissue decay, an asymptotic analysis of (18)–(20) for this case is presented in Sect. 5. Overall, the results in this section suggest that must hold for tissue growth to occur, thus suggesting that tissue engineers should ensure that the porosity of the scaffold is at least larger than the ratio between the rate of cell death and growth.

Travelling-Wave Analysis for

Figure 1a, b indicates the emergence of semi-infinite travelling waves of constant speed for In light of this, we assume that for sufficiently large time, where c is the constant wave speed at which the tissue edge moves. In this section, we employ travelling-wave analysis to obtain the wave speed c in terms of the governing parameters when .

We transform (18)–(20) via the travelling-wave coordinates where . Setting , we obtainFollowing Fadai and Simpson (2020) and Fadai (2021) , we defineMultiplying (25) by and rewriting the conditions from (26) in terms of q(z), we obtainHere, we note the second boundary condition from (26) transforms into the second of (29) because Dividing (28) by (27) we haveUsing the shooting method (presented in “Appendix B”) to find the heteroclinic connection q(n) that connects to , we can determine a numerical approximation of the wave speed in terms of the governing parameters and

In Fig. 2a, b, the solid black line represents the relationships and c(s), respectively, when (30) and (31) is approximated by the shooting method. The dashed green line represents these wave speeds when obtained by numerically solving (18)–(20), and computing c by evaluating at large time. In view of the close agreement between these two approaches to computing c,  we henceforth concentrate on solutions provided by (30) and (31) for simplicity.

Fig. 2

Numerical solutions for the wave speeds and c(s). The solid black and dashed green lines represent approximations sought by numerically solving (30) and (31) using a shooting method and the system from (18)–(20), respectively. Parameter values for a and Parameter values for b and (Color figure online)

The results presented in Fig. 2 suggest that larger and smaller s increase the speed at which the tissue front grows, and since this further corresponds to minimising and maximising and the porosity of the scaffold. In Fig. 3a, b we present the wave speeds for and respectively. These results suggest that, for a fixed value of , the wave speed is maximised when Physically, this corresponds to the case where the viscosity of the cells is much greater than the viscosity of the liquid, and where repulsive forces exerted due to cell–cell and cell–scaffold interactions at high cell volume fractions dominate inter-cellular forces that give rise to cell dispersal. Furthermore, Fig. 3a, b indicates that the dependence c has on is weaker for than This suggests that for smaller , cell–cell and cell–scaffold interactions which can drive cellular motion are more prominent in controlling the wave speed than the cell and liquid viscosities. Additionally, and in agreement with Fig. 2a, Fig. 3 indicates that the wave speed increases as increases.

Fig. 3

Numerical solutions for the wave speeds for and sought by numerically solving (30) and (31) using a shooting method. was used for both sub-figures. We note that the colour axes are different in each sub-figure (Color figure online)

Asymptotic Analysis for

In this section, we construct asymptotic solutions for and L(t) for when Since scaffold porosity is a readily controllable parameter (in contrast to cell growth and death), the analysis in this section when can be associated with the case in which the scaffold porosity is low and tissue growth is successful, so that

The numerical results in Fig. 1e indicate that holds for when Furthermore, given that when and is also expected when for sufficiently large time. We therefore have when and . To aid the subsequent asymptotic analysis, we introduce the variableswhich imply that (18) and (19) are simplified toTo analyse the behaviour of (33) and (34), we follow Newman (1980) and adopt the ansatzwherein Imposing on (35), we obtain thatIn view of (35), (33) providesIn this section, initial conditions for N and L are chosen to satisfy (35) and (36) when – i.e., and hence An implicit solution for N and L is found by computing A(B) from (37) and exploiting the relation from (36). Following Newman (1980), we havewhere and

We now deduce the large-T behaviour of and L(T),  from which the large-time behaviour of and L(t) when can subsequently be determined. Guided by the numerical results from Fig. 1d, the evolution of L satisfies for sufficiently large T, so we have from (39) thatwhich is then inverted to givefor . Equation (38) and the leading-order term in the above expansion are used to find the following large-T approximation for N:We now exploit (42) and (41) to deduce the large-time behaviour of and L(t) when

Large-Time Behaviour of n and L When

When , then takes the form of a travelling wave of constant speed and for In contrast, when , the numerical results in Fig. 1f suggest that as for some finite constant In general, a large-time solution for L is unavailable given that does not necessarily hold when however, when then to leading order at and L evolves according to (41) until . Since when and , (41) becomesfor and Equation (43) implies that when and Therefore, travelling waves propagate with speed when this being in agreement with the numerical results from Fig. 2a. For , (43) implies that the growth of the tissue edge is logarithmic until and linear thereafter. If tissue growth is successful, this suggests the formation of travelling waves with constant speed is delayed when the scaffold porosity is low. When , the exponential terms from (43) are negligible as and we obtainwhen and The leading order logarithmic term in (43) and (38) is used to find the following large-time approximations for n when We note that (43) and (45) hold for if and because and in this case. Furthermore, for and then (38) suggests that at large time since

In Figs. 4 and 5, we compare the numerical solution for n and L when obtained by numerically solving the PDE system from (18)–(20) for and against their respective asymptotic solutions from (45), (46) and (43). Overall, an excellent agreement between the numerical and asymptotic solutions is observed. The large-T behaviors of N and L characterised by (41) and (42), and hence the asymptotic approximations from this subsection, are only valid for initial conditions that satisfy (35) and (36). We now show that solutions of (33) and (34) converge to solutions similar to that of (41) and (42) for a wider class of initial data.

Fig. 4

Numerical solution (solid black line) of the PDE system from (18)–(20) versus asymptotic solution (dashed green line) from (45) and (43) for n (a) and L (b) for . Solutions for n are presented on at intervals of 1000. Initial conditions for the numerical simulations were chosen to satisfy (35) and (36). Parameter values: and (Color figure online)

Fig. 5

Numerical solution (solid black line) of the PDE system from (18)–(20) versus asymptotic solution (dashed green line) from (46) and (43) for n (a) and L (b) for . Solutions for n are presented on at intervals of 500. Initial conditions for the numerical simulations were chosen to satisfy (35) and (36). Parameter values: and (Color figure online)

Convergence of Asymptotic Solutions

Since the choice of initial cell distribution within the scaffold is likely to vary substantially in practice, it is important to determine the large-T behaviour of N and L for a wider class of initial data, such as those from (22). The asymptotic behaviour of (33) as comprises an interior layer near the interface within which the similarity reduction holds where and so that satisfiesFollowing familiar arguments to that of the porous-Fisher equation (Aaronson 1980; Murray 2002), although we emphasise that N is not a travelling wave of the usual form, the solution to (47) is given byFor settingimplies thatthe solution to which that matches into the exponential terms in (48), and hence the corresponding term in the interior layer, is given byfor some unknown constant m, and this dominates the asymptotic behaviour of (50) as . Therefore, (49) becomeswhere m depends on the initial data. Since (52) implies thatfor . By comparing (52) and (42), we see that the asymptotic structure is retained despite the initial cell distribution for large T. In addition, (53) suggests that the choice of initial cell distribution does not affect the speed at which the tissue edge moves for large T, but does affect the position of the tissue boundary. We note that if is chosen to satisfy (35), then (42) indicates that

In Fig. 6, we compare the numerical solutions for N and L when obtain by numerically solving the PDE system from (33) and (34) against the asymptotic solutions from (52) and (53) for two choices of and . The value of was found by solving For both , we are able to choose an m that provides excellent agreement between the numerical and asymptotic solutions. We note that the large-time behaviour for and L(t) when can be extracted directly from (52) and (53) given that when and This justifies the numerical results observed in Fig. 1c, d.

Fig. 6

Numerical solution (solid black line) of the PDE system from (33) and (34) versus asymptotic solution (dashed red line) from (52) and (53) for n and L for different initial conditions and . Solutions for n are presented on at intervals of 1000 (Color figure online)

Conclusions

In this paper, we present a multiphase model to represent the development of tissue in vitro within a porous tissue engineering scaffold. We consider a tissue cell phase, extra-cellular liquid phase and a scaffold phase and adopt Darcy’s law to relate the velocity of the cell and liquid phases to their respective pressures. The model includes mechanisms to represent cell growth and death, and pressures that arise from cell–cell and cell–scaffold interactions. We employ a moving boundary, to track the speed at which the tissue edge propagates through the scaffold. We reduce the model to a nonlinear reaction–diffusion equation for the cell volume fraction, , and a moving boundary condition for the tissue edge. The diffusivity of the reaction–diffusion equation is dependent on the cell and scaffold volume fractions; cell and liquid viscosities, and pressures that arise from cell–cell and cell–scaffold interactions. Non-dimensionalisation of the model shows that the tissue boundary position increases with the scaffold permeability and exposes important dimensionless groupings. One such grouping, that describes the difference between the scaffold porosity and the ratio between the cell death and growth rates is of crucial importance to the qualitative features of the cell phase evolution. The model admits three regimes for the evolution of the cell volume fraction and the moving boundary, based on the sign of . Employing travelling-wave and asymptotic analysis, we characterise these regimes in terms of and parameters related to cellular motion.

The case in which corresponds to the successful growth of tissue, which suggests that tissue engineers should ensure that the porosity of the scaffold is at least larger than the ratio between the rate of cell death and growth. For we show that the cell volume fraction, spreads through the scaffold as a semi-infinite travelling wave with constant speed, emerging from the steady state Employing travelling-wave analysis, we accurately compute the wave speed (i.e., the speed at which the tissue edge moves through the scaffold) as a function of the governing parameters. We find that the wave speed is greatest when the rate of apoptosis is negligible in comparison with that of mitosis, the viscosity of the cells is much greater than the viscosity of the liquid, and when repulsive forces exerted due to cell–cell and cell–scaffold interactions at high cell volume fractions dominate inter-cellular forces that give rise to cell dispersal. We also find that the wave speed increases as the scaffold porosity increases; however, we note that the cells will require a sufficient amount of scaffold on which to attach, so an upper bound on the porosity is to be expected. Furthermore, we deduce that for smaller values of and hence scaffolds with small porosities, cell–cell and cell–scaffold interactions which can drive cellular motion are more prominent in controlling the wave speed than the cell and liquid viscosities.

For we employ asymptotic analysis to find explicit solutions for n and L. Since scaffold porosity is a readily controllable parameter (in contrast to cell growth and death), the analysis in this section when can be associated with the case in which the scaffold porosity is low and tissue growth is successful, so that When the growth of the tissue edge is logarithmic until and linear thereafter, thus suggesting the formation of travelling waves with constant speed is delayed as and hence when the scaffold porosity is low. For we deduce that the cell volume fraction decays exponentially with rate at large time, with the moving boundary tending towards a steady state. For and , the evolution of the L is shown to be logarithmic until and approaches a steady state thereafter, the value of which is found explicitly and related to and the initial conditions employed in the model. For we also demonstrated that the choice of initial cell distribution does not affect the eventual distribution of cells within the scaffold, nor the speed at which the tissue edge moves, but does affect the position of the tissue boundary.

For a functional tissue construct to develop within a scaffold, cells must be exposed to the correct environment and stimulated with growth factors such as oxygen. There must also be a sufficient amount of scaffold on which the cells can adhere. Whilst key features of tissue growth such as cell mitosis, apoptosis and motion are included in this paper, concepts such as environmental pressures, cellular adhesion, and nutrient supply have not been considered. Therefore, following Lemon and King (2007a), a natural extension of this work would include examining the influence that nutrient limitation has on cell growth. We leave these extensions for future consideration.