Numerical Methods for Two-Dimensional Stem Cell Tissue Growth.

Growth of developing and regenerative biological tissues of different cell types is usually driven by stem cells and their local environment. Here, we present a computational framework for continuum tissue growth models consisting of stem cells, cell lineages, and diffusive molecules that regulate proliferation and differentiation through feedback. To deal with the moving boundaries of the models in both open geometries and closed geometries (through polar coordinates) in two dimensions, we transform the dynamic domains and governing equations to fixed domains, followed by solving for the transformation functions to track the interface explicitly. Clustering grid points in local regions for better efficiency and accuracy can be achieved by appropriate choices of the transformation. The equations resulting from the incompressibility of the tissue is approximated by high-order finite difference schemes and is solved using the multigrid algorithms. The numerical tests demonstrate an overall spatiotemporal second-order accuracy of the methods and their capability in capturing large deformations of the tissue boundaries. The methods are applied to two biological systems: stratified epithelia for studying the effects of two different types of stem cell niches and the scaling of a morphogen gradient with the size of the Drosophila imaginal wing disc during growth. Direct simulations of both systems suggest that that the computational framework is robust and accurate, and it can incorporate various biological processes critical to stem cell dynamics and tissue growth.

Introduction

During development and regeneration of biological tissues, stem cells play a crucial role in controlling growth, morphology, and tissue size. Stem cells proliferate to maintain their population while they differentiate to give rise to other types of cells for specification of their fates, leading to tissues of different biological functions and properties. The stem cell population and tissue growth are intimately connected, and stem cell behavior is governed by many intrinsic and extrinsic factors, often through regulation of the self-renewal and the cell cycle length of stem cells [33]. Such regulatory factors include intracellular signaling, diffusive molecules secreted by cells [15], cell-to-cell interactions [4], mechanical forces, and contact-mediated intercellular signaling of morphogens [15]. As a result, the stem cell niche, the microenvironment in which stem cells reside in a tissue, has been identified as a key factor for the sustenance of stem cells and regenerative tissues as a whole [19, 37, 43].

To model spatial and temporal dynamics of stem cells and tissue growth, the discrete cell approach, such as cellular automata [57] or subcellular element methods [44], may be used directly to track each individual cell in the tissue for its growth and spatial dynamics. Such techniques are computationally convenient and useful when the number of cells is small or when the individuality of stem cells and their discrete nature are important in modeling cell lineages and growth [12, 49]. For tissues of large numbers of cells, continuum models, in which cell densities are functions of space and time, are governed by conservation laws, such as those for the spreading and invasion of solid tumors [2, 13, 24, 25], may be more appropriate. This avenue provides a robust framework that allows for direct incorporation of many important biochemical regulations and physical and mechanical components often observed in stem cell niches and tissue growth.

Recently, a one-dimensional continuum model has been developed to study the spatial aspects of cell lineages and the formation and sustenance of the stem cell niche in the olfactory epithelium (OE) [11]. Using morphogens to provide feedback on the fates of stem cells together with a permeable basal lamina, the model achieves a stratification of cell layers as the tissue grows to a homeostatic state during development. In this model, the geometry of the tissue is simplified to one spatial dimension along the basal-apical axis with one end of the tissue at the basal lamina assumed to be fixed while the other end is dynamic due to proliferation and differentiation of stem cells and other cell types in the tissue. While the model is focused on the stratification of different types of cells in the vertical direction, it neglects influences from neighboring cells in the horizontal direction parallel to the basal lamina that have been observed experimentally to exhibit strong spatial inhomogeneity in many types of epithelia [20, 22, 35].

Continuum tissue modeling in high dimensions, as opposed to one dimension, presents significant computational challenges. Often some form of a mechanical closure is required to solve for the tissue growth velocity. Because tissue growth involves interfacial motion of boundaries of the tissue domain, the two-dimensional models also require special computational techniques for the treatment of such moving boundaries. Level set methods [47, 53], phase field methods [9, 61], and immersed interface methods [31, 39] are robust approaches, particularly suitable for dealing with topological changes of the dynamic interfaces. Interface tracking methods, such as the boundary integral approach [38] and front-tracking methods [21, 40, 56], are efficient and accurate in capturing fine structures of the interfaces.

For tissue growth driven by a stem cell niche, several crucial aspects must be taken into consideration when choosing numerical methods for interfacial motion. Often, intertissue signaling of diffusive regulatory molecules occurs across tissue boundaries [11], indicating that methods capturing the boundary in a non-diffuse fashion may be most appropriate. The presence of cell-to-cell adhesion [24] may intricate details of the tissue boundary and require accurate tracking of its movement. To capture long-term distortions of tissue morphology, preservation of the interface over long time scales may be necessary. Also, accounting for topological changes of the tissue domain may not be of particular importance as budding, pinching, and fusion of domains rarely occur with developing and regenerative tissues.

Here, we introduce numerical methods for stem cells, cell lineages and tissue growth in two spatial dimensions through a continuum modeling framework. We develop an interface-tracking method based on a transformation technique recently developed for studing viscous effects of Stokes waves using the incompressible Navier–Stokes equations in two spatial dimensions [58]. In this approach, a single interfacial boundary between two fluids is dealt with by transforming the deformed geometry and governing equations to the unit square. The resultant incompressibility equation with variable coefficients are then solved by an iterative procedure that enables usage of the pseudospectral method in the horizontal direction in which periodicity is assumed. This method, which explicitly tracks the interface, can easily incorporate transformations that allow clustering of points in localized spatial regions for better accuracy and efficiency.

Our models for growing tissues of different cell types consist of convection equations for cell movement and proliferation, equations for tissue incompressibility and Darcy’s law, and reaction–convection–diffusion equations for diffusive molecules secreted by cells. At the center of our numerical approach is a high-order finite difference approximation for solving the tissue incompressibility equation along with a multigrid solver. To handle more complex growth features, we extend the one dynamic boundary transformation-based approach to domains with two dynamic boundaries. We also re-formulate the model in polar coordinates to handle closed geometries and introduce a corresponding transformation to the unit circle for the usage of multigrid methods in polar coordinates. The convection terms are treated with the second-order upwind schemes and the temporal discretization is approximated by second order TVD Runge–Kutta methods [23].

Numerical tests of various cases demonstrate that our computational framework is second-order accurate in both time and space with one or two dynamic boundaries in rectangular coordinates and in polar coordinates. We also demonstrate through examples that the combination of high-order finite difference spatial discretizations in conjunction with multigrid solvers provide a robust method for the movement of tissue interfaces after transformation. In particular, it is shown that this approach is more efficient and robust in capturing significantly deformed interfaces than the pseudospectral approach with a direct iterative solver [58].

To apply our methods to biological tissues, we use the stem cell niche in stratified epithelia, an open geometry system, and the Drosophila imaginal wing disc, a closed geometry system, as models. We first extend the one-dimensional OE model [11] to general stratified epithelia with two different types of stem cell niches in two dimensions (this model will be further studied in [48]). The imaginal wing disc serves as an application for our methods in polar coordinates to study the scaling of morphogen gradients during growth.

The remainder of this paper is organized as follows. In Sect. 2, we construct methods to solve equations governing tissue growth on a moving domain with one and two dynamic boundaries in rectangular coordinates; in Sect. 3, we develop such methods for polar coordinates; in Sect. 4, we provide computational results to test our methods and compare them to an iterative pseudospectral approach; in Sect. 5, we apply our methods to the stem cell niche in stratified epithelia and to the Drosophila imaginal wing disc; and in Sect. 6, we conclude the paper.

Numerical Methods in Rectangular Coordinates

One Dynamic Boundary

Consider the two-dimensional domain as a spatial representation of a biological tissue. We allow and to represent a dynamic free-top boundary and fixed bottom of the tissue, respectively. These two boundaries may be considered to be either the interface between neighboring tissues that may provide external input into the system or an open external environment. We will assume periodicity in the -direction in throughout, essentially modeling just one segment of a wide, encapsulating tissue. We also take to be the velocity of the tissue.

Suppose that the tissue is composed of different types of cells. Assuming that the diffusion of cells is negligible, the equations that govern the densities of each cell type, for , are,where is the net rate of production (or removal) of the th cell type. Further assuming a uniform density of cells within the tissue, or that the tissue behaves as an incompressible fluid, we have , and, therefore,To solve for , we introduce a new variable representing pressure. If we assume that the tissue behaves similar to a porous medium, then we may relate pressure and velocity by Darcy’s law [6, 60],which generally states that tissue flow occurs down the pressure gradients. Combining Eqs. (2) and (3) yieldsTo maintain as a fixed, immobile boundary, we require the vertical velocity component to be zero there,Along the boundary at the top of the tissue, , we assume there exists a surface tension that may result from intercellular forces and cell surface mechanics [17, 36]. Thus, we take the pressure to be directly proportional to the curvature of the boundary, , [25],The size of the proportionality constant corresponds to the magnitude of the surface tension that exists along , which may vary depending on the number and strength of adherens and gap junctions between cells. The growth of the dynamic boundary of the system, , is governed by the kinematic condition,which requires the normal component of the tissue velocity at the boundary to be equal to the velocity of the tissue boundary itself.

Transformation to the Unit Square

To handle the governing partial differential equations, we first transform the domain to the unit square and then track the movement of the dynamic boundary explicitly. Following [58], we apply the following coordinate system, where and . The transformation function can be chosen to naturally accumulate grid points near the dynamic interface where morphogen gradients are sharp. The transformed derivatives are written as follows: whereFurthermore, we can write the Laplacian operator aswhere

Scaling the Tissue Growth System

For simplicity, we first scale the by ,Applying the transformation to Eqs. (3–4) yields The boundary conditions in Eqs. (5–6) reduce to where . The description of the movement of the boundary in Eq. (7) is essentially invariant under the transformation,

Two Dynamic Boundaries

Let us also consider the case in which the tissue does not have a rigid boundary at , which is instead replaced by a second dynamic boundary. In this case, we may begin with the domain to represent a biological tissue where both and are dynamic boundaries. If we keep to our previous assumptions of cells comprising the tissue with a uniform density and impose Darcy’s law, then the governing equations are still given by Eqs. (1–4). Continuing to assume that surface tension is is present along the dynamic boundaries of the domain, we have the following boundary conditions, where and are the respective curvatures of and . Also, the kinematic boundary conditions are given by, Transforming the system on to a rectangular domain occurs in a similar fashion as described in Sect. 2.1.2 by choosing and .

Spatial and Temporal Discretizations and Multigrid Solvers

Now that the governing equations on and have been transformed to a rectangular domain, we then present methods to solve the transformed system in Eqs. (20–24) for the case with one dynamic boundary and for the case with two boundaries presented in Sect. 2.2.

Both the curvature and are functions of in the transformed domain regardless of the particular choice of in Eq. (9), since takes two spatial derivatives of . As a result, depends on while the partial derivatives of govern the movement of as evident in Eq. (7). Naturally, it follows that the movement of is implicitly dependent upon its own third derivative, , and solving the transformed system of equations requires a high-order accuracy of the spatial discretization in the -direction to maintain overall second-order accuracy in at least computing . Initial numerical tests in which the system was discretized in space using second-order central difference approximations and then evolved in time were observed to be unstable and inaccurate.

One approach is to use pseudospectral method in the -direction similar to methods for the transformed inviscid Navier–Stokes equations in conjunction with an iterative method [58]. In Appendix 7, we outline an algorithm for such an approach to solve Eqs. (20–24) for a single dynamic boundary. However, an iterative algorithm incorporating the pseudospectral method has not yet been applied to domains with two dynamic boundaries. We also note that the required numerical accuracy for these equations are not as high as those for the inviscid Navier–Stokes equations. In Sect. 4, we will present the performance of such an implementation.

Here, we use fourth-order central difference approximations for the derivatives of and the partial derivatives of in the -direction, namely, Because of the assumed periodic boundary conditions in the direction of and in the system, the five-point stencil in the -direction near the boundary can be implemented in a straightforward fashion. For the discretization of the partial derivatives in the -direction, we use second-order central difference approximations.

To solve the linear system that results from the spatial discretization of Poisson’s equation in Eq. (21), we implement a Gauss–Seidel iterative scheme with a multigrid solver, though other schemes, such as Newton relaxation, may be used. Our transformed Laplacian operator described in Eq. (17) stands very different from the usual Laplacian operator on the unit square due to first-order and mixed partial derivatives that are present, and, more importantly, variable coefficients for three of the terms. This notion draws similarities between our transformed Laplacian operator and the anisotropic diffusion operator. Several multigrid algorithms have been constructed to robustly and efficiently handle the anisotropic diffusion operator and general cases with variable coefficients [1, 16, 42, 46, 51], which may be applied here. However, for simplicity, we implement a multigrid algorithm that is often applied to solve elliptic problems with constant coefficients discussed in [8] with a 9-point full-weighting for the restriction operator. As a result of the variable coefficients, the number of Gauss–Seidel relaxations we would need to implement on each level to solve the linear system exactly is dependent upon the size of the original grid, which stands in contrast to the case with constant coefficients. In all of our numerical simulations, we use 80 Gauss–Seidel relaxations on each level to ensure the multigrid algorithm solves the resultant linear system accurately, though this many relaxations on each grid is often unnecessary.

To explicity track the movement of the boundary governed by the kinematic equation in Eq. (24), or for and in the case of two dynamic boundaries, a second-order upwind approximation is used to discretize in spaceFor time evolution, a second-order TVD Runge–Kutta is employed [23], with corresponding multigrid solvers used during each corresponding function evaluation, . As a result, the overall spatiotemporal order of accuracy of this computational approach is second-order.

Numerical Methods in Polar Coordinates

In this section, we translate the moving boundary system in rectangular coordinates described in Sect. 2 to a formulation in polar coordinates. We begin with a two-dimensional polar domain to represent our tissue where is a dynamic free-boundary. We take to be the tissue velocity. If we assume that cell types compose the tissue as an incompressible fluid that obeys Darcy’s law, then we again have the governing equations, where variables are defined similarly as in our rectangular formulation. The kinematic equation describing the dynamics of is

Transformation to the Unit Circle

We transform the system in Eqs. (34–37) from to the unit circle similar to the manner in which is scaled to the unit square: where and . The transformed derivatives are written as follows: withThe transformed Laplacian takes the form,where

Scaling the Tissue Growth System in Polar Coordinates

For the governing equations, we first scale by ,Applying the transformation to Eqs. (34–35) yields The boundary condition in Eq. (36) reduces towhere . The kinematic condition in Eq. (37) scales as

Spatial and Temporal Discretization and Multigrid Solvers in Polar Coordinates

In polar coordinates, the overall computational approach stands very similar to the multigrid approach for rectangular coordinates developed in Sect. 2. A fourth-order central difference approximation is used to discretize in the -direction due to the dependence of and on while second-order finite difference approximations are used in the -direction. The resultant linear system is then solved by a Gauss–Seidel iterative scheme with a multigrid solver. For interpolation on the discretized unit circle, we use inverse distance weighting that depends on . More specifically, our interpolation operator from takes the form, The restriction operator used is just the transpose of the interpolation operator scaled by the normalization constant , which results in an 11-point full-weighting. Due to variable coefficients in the transformed Laplacian operator in polar coordinates, we again use 80 Gauss–Seidel relaxations on each grid for the general multigrid algorithm, though multigrid methods specialized for polar coordinates [5, 28, 32] may be applied here.

The behavior of the system at the singularity at the origin depends on the choice of that defines our transformation for the polar coordinate system. Choosing a linear scaling of results in all the coefficients of the transformed Laplacian being directly proportional to inverse powers of ,This form, then, allows for a simple treatment of the Laplacian operator at the origin,where is the spatial discretization size in the -direction.

For movement of the dynamic boundary in polar coordinates and temporal discretization, upwind discretizations and Runge–Kutta methods are used as described for rectangular coordinates.

Numerical Tests

In this section, we demonstrate that the order of accuracy of our method for interfacial motion using a transformation achieves an overall second-order spatiotemporal accuracy. We also compare the performance of implementing a multigrid solver with a fourth-order central difference approximation in the -direction with the iterative pseudospectral approach outlined in Appendix 7.

To compare these two methods, we implement both of them to solve an exact static problem for the pressure, , without time evolution given by Eqs. (4–6). The variable boundary along the top of the tissue is described by with varying in size for the numerical tests. We take a simple linear form of in Eq. (9) for our transformation in this testing case, implying grid points are evenly distributed in the -direction. For , we see that the iterative pseudospectral approach achieves a second-order accuracy with respect to the maximum-norm and with approximately 10–20 iterations required for convergence as evident in Table 1. However, Table 1 also shows that when the value of is increased to 0.15, many more iterations are required for convergence on grids of size and 16, while the method fails to converge on the larger grids tested. The iterative pseudospectral approach fails to converge in our numerical simulations on any of the grid sizes tested when . We speculate that this is due to the increased magnitude and variability of the coefficients of the transformed Laplacian operator in Eq. (17).

Table 1

Errors, orders of accuracy, and CPU times for calculation of the internal tissue pressure for the case with one variable boundary given by Eqs. (4–6) without time evolution using the iterative pseudospectral approach in Appendix 7

N	B=0.1				B=0.15				B=0.25
	Error	Order	Iter.	CPU(s)	Error	Order	Iter.	CPU(s)	Error	CPU(s)
8	0.466	–	16	5.40e−2	2.17	–	34	9.90e−2	NC	–
16	0.100	2.22	19	0.182	0.226	3.27	64	0.226	NC	–
32	2.55e−2	1.98	21	0.624	NC	–	–	–	NC	–
64	6.42e-3	1.99	21	2.19	NC	–	–	–	NC	–
128	1.61e−3	2.00	18	6.93	NC	–	–	–	NC	–
256	4.02e−4	2.00	15	20.4	NC	–	–	–	NC	–
512	1.00e−4	2.00	12	63.9	NC	–	–	–	NC	–

The boundary is described by the curve . An exact testing case is constructed with and a corresponding where is the curvature of . A linear scaling is used for the transformation in Eqs. (8–10) along with a tolerance of for the iterative method. Parameters chosen are

NC no convergence of the iterative method

Table 2 demonstrates that a multigrid solver with a fourth-order central difference approximation in the -direction also approximates with a second-order spatial accuracy in the case when . Comparing the corresponding CPU times along with those for the iterative pseudospectral simulations reveals that, even with 80 Gauss–Seidel relaxations on each grid, the multigrid approach is the less expensive of the two algorithms. Namely, the multigrid tests achieve a twofold increase in CPU time for the largest grid tested with . As increases to and in our tests, Table 2 displays convergence and a second-order accuracy for the multigrid approach with no increase in CPU time due to the fixed number of Gauss–Seidel relaxations. This notion stands as a drastic improvement over the pseudospectral tests in which the iterative solver often failed to converge when was further distorted with values exceeding 0.1. We note that some accuracies computed for this approach achieve orders of accuracies exceeding two with a highly deformed interface, and we speculate that this occurs due to a dominance of the error in the spatial direction in which fourth-order central difference approximations are made.

Table 2

Errors, orders of accuracy, and CPU times for calculation of the internal tissue pressure for the case with one variable boundary given by Eqs. (4–6) without time evolution using fourth-order central difference approximations along with a multigrid solver

N	B=0.1			B=0.25			B=0.5
	Error	Order	CPU(s)	Error	Order	CPU(s)	Error	Order	CPU(s)
8	0.512	–	1.30e−2	29.0	–	1.50e−2	412.8	–	1.20e−2
16	6.23e−2	3.04	1.70e−2	1.99	3.87	1.70e−2	128.2	1.69	1.50e−2
32	1.22e−2	2.35	6.60e−2	7.54e−2	4.72	6.50e−2	5.86	4.45	6.40e−2
64	2.78e−3	2.14	0.256	9.65e−3	2.97	0.256	0.168	5.61	0.242
128	6.69e−4	2.06	1.12	1.69e−3	2.51	1.15	8.97e−3	3.75	1.12
256	1.64e−4	2.02	5.92	4.19e−4	2.01	5.84	1.04e−3	3.11	5.99
512	4.08e−5	2.01	29.5	1.05e−4	2.00	30.6	2.74e−4	1.92	30.2

The variable boundary , choice of transformation, and parameters chosen are the same as those given in Table 1

From here, we proceed to test the accuracy of the overall spatial and temporal discretizations. Table 3 shows an overall second-order spatiotemporal accuracy for the coupled system of both and using maximum-norm difference between successive approximations to compute the error and, as a result, to compute the order of accuracy. Since moves dynamically in this testing case, an exponential scaling of in Eq. (9) is used for a higher spatial resolution in the -direction near . The transformed grid with is plotted on the initial state of in Fig. 1a. We also demonstrate a proper temporal treatment of the kinematic boundary condition in Eq. (24) in Table 4 using an exact dynamic problem in which is spatially uniform.

Table 3

Errors, orders of accuracy, and CPU times for calculation of the internal tissue pressure and boundary movement for the case with one dynamic boundary given by Eqs. (4–6) up to time

N	P		h		CPU(s)
	Error	Order	Error	Order
8	–	–	–	–	3.60e−2
16	4.28e−2	–	8.44e−3	–	0.504
32	1.022e−2	2.07	2.41e−3	1.81	8.13
64	2.54e−3	2.01	5.35e−4	2.17	136.0
128	6.27e−4	2.02	1.23e−4	2.12	2,372.3
256	1.55e−4	2.02	2.97e−5	2.05	46,870.3

The initial state of dynamic boundary is described by the curve , and a uniform influx of cells is assumed by . An exponential scaling is used for the transformation in Eqs. (8–10) with . is taken for time discretization. Parameters chosen are and

Fig. 1

plotted on computational grids for numerical tests in rectangular coordinates with shown on the domains and . a An exponential scaling for the transformation in Eqs. (8–10) with and as implemented for calculations in Table 3 and Fig. 3. b An arctangent scaling for the transformation in Eqs. (8–10) with and and as implemented for calculations in Table 5

Table 4

Errors, orders of accuracy, and CPU times for calculation testing the implementation of the kinematic boundary condition for the case with one dynamic boundary given by Eqs. (4–6) up to time

Δτ	P		h		CPU(s)
	Error	Order	Error	Order
2.0e−2	1.71e−2	–	1.89e−4	–	4.69
1.0e−2	4.46e−3	1.94	4.91e−5	1.94	9.28
5.0e−3	1.14e−3	1.97	1.26e−5	1.96	18.7
2.5e−3	2.87e−4	1.99	3.16e−6	2.00	37.7
1.25e−3	7.20e−5	2.00	8.33e−7	1.92	75.1

The initial state of the dynamic boundary is described by , and a uniform influx of cells is assumed by . The exact solutions are given by and . A linear scaling is used for the transformation in Eqs. (8–10). is taken for spatial discretization. Parameters chosen are and

Fig. 3

Distributions of the densities of stem cells (), TA cells (), TD cells (), concentrations of and , and pressure () at time for stratified epithelium simulations with a free-form stem cell niche. Initial conditions and the assumed spatiotemporal discretization are those used for this figure, and parameters used are given in Table 7

Table 5

Errors, orders of accuracy, and CPU times for calculation of the internal tissue pressure and boundary movement for the case with two dynamic boundary given by Eqs. (4, 25–28) up to time

N	P		h_1		h_2		CPU(s)
	Error	Order	Error	Order	Error	Order
8	–	–	–	–	–	–	3.10e−2
16	9.42e−3	–	7.53e−3	–	2.04e−2	–	0.427
32	3.71e−3	1.34	9.47e−4	2.99	6.65e−3	1.61	7.83
64	1.11e−3	1.74	1.44e−4	2.72	2.40e−3	1.47	123.6
128	2.75e−4	2.02	3.10e−5	2.21	5.67e−4	2.08	2,344.7
256	6.96e−5	1.98	8.55e−6	1.86	1.46e−4	1.95	47,573.1

The initial states of the dynamic boundaries and are described by the curves and , and a uniform influx of cells is assumed by . An arctangent scaling is used for the transformation in Eqs. (8–10) with . is taken for time discretization. Parameters chosen are and

We next test the accuracy of our approach for the problem formulated with two dynamic boundaries in Sect. 2.2. Table 5 shows that approximations for and achieve a second-order accuracy in a numerical testing case evolved in time using the maximum-norm difference between successive approximations to compute the error. For this test using two variable boundaries, an arctangent scaling of is used to place more grid points in the -direction near both and . The transformed grid is plotted on the initial state of for the dynamic problem in Fig. 1b.

To test the accuracy of our approach with a fourth-order central difference approximation in the -direction in polar coordinates, we choose a dynamic time-dependent problem with no exact solution and a constant form of . A second-order accuracy in time and space is achieved using differences between successive approximations to compute the error as Table 6 demonstrates. As previously discussed, a linear form of in Eq. (38) is chosen for ease in handling the singularity at the origin.

Table 6

Errors, orders of accuracy, and CPU times for calculation of the internal tissue pressure and boundary movement in polar coordinates given by Eqs. (4–6) up to time

N	P		H		CPU(s)
	Error	Order	Error	Order
16	–	–	–	–	0.152
32	2.47e−3	–	5.90e−3	–	2.70
64	3.97e−4	2.63	1.50e−3	1.98	45.5
128	9.38e−5	2.08	3.85e−4	1.96	895.9
256	1.96e−5	2.26	8.97e−5	2.10	26,032.1

The initial state of dynamic boundary is described by the curve , and a uniform influx of cells is assumed by . A linear scaling is used for the transformation in Eqs. (38–40). is taken for time discretization. Parameters chosen are and

Applications

In this section, we apply our numerical methods for tissue growth to two model systems in developmental biology: stratified epithelia and the Drosophila imaginal wing disc. The application to developing stratified epithelia focuses on the spatial aspects of cell lineages and stem cell niche formation via gradients of signaling molecules in open geometries using the two-dimensional rectangular coordinate system approach presented in Sect. 2. For tissue growth in the closed geometry of the imaginal wing disc, we test our methods in polar coordinates for an expansion–repression mechanism on a growing domain and its effects on the scaling of morphogen gradients. For both systems, proliferation and its interplay with biochemical signaling processes play critical roles in tissue growth.

A Spatial Cell Lineage Model in Two-Dimensions

Stratified epithelia are great model systems to study the development, sustenance, and regeneration of tissues through a cell lineage, in which the progeny of cells progress from stem to terminal differential stages. Many types of epithelia that are several cells thick stratify into different cell layers and maintain a population of stem cells along the basal lamina of the tissue. Examples of such tissues are evident in the epidermis [30], OE [54], and cerebral cortex [18]. In these tissues, the stem cell niche plays the crucial role of controlling the populations of cell types that stratify in layers during developmental, homeostatic, and regenerative conditions [50]. Though the patterning of cell layers occurs primarily along the apical-basal direction, spatially varying two- and three-dimensional tissue morphologies have been observed in stratified epithelia [20, 22, 35]. Here, we will extend the one-dimensional mathematical model for the stratification of cell layers and the formation of a stem cell niche in the OE [11] to a model for stratified epithelia in two-dimensions [48].

Assume that the epithelium is represented by the domain , and it is comprised by three different cell types: stem cells, transit amplifying (TA) cells, and terminally differentiated (TD) cells. Representing the densities of each of the cell types by , and , respectively, the governing equations are where is the replication probability of the th cell type and is over its cell cycle length for . The natural death rate of TD cells is represented by . It then follows that the net influx of cells into the tissue is given byWe also assume no-flux boundary conditions for each of the cell densities along the top and bottom of the tissue. To induce stratification of the tissue into layers, we allow morphogens that are produced by cells in the epithelium to be rapidly uptaken into the underlying basement membrane to form spatial gradients and influence the probabilities and . If we assume that the diffusive molecules and are produced by cells at a rate of and , respectively, then both morphogens are described by Here, and are diffusion coefficients and and are degradation rates for each of the molecules. To interpret the rapid uptake of each of the molecules, we employ a leaky boundary condition for each along the basal lamina and a no-flux boundary at the closed apical surface. If we take the bottom of the tissue, , to be the basal lamina and the top as the apical tissue surface, then the boundary conditions become where denotes the unit outward vector normal to and and represent the coefficients of permeability of each of the respective morphogens as they diffuse through the basal lamina. The stem cell niche that develops in this case is referred to as the rigid stem cell niche since it prompts niche formation along the rigid boundary . If the top of the tissue, , is taken as the basal lamina and the bottom as the apical surface, then our boundary conditions take the form, This case is referred to as the free-form stem cell niche due to niche formation along the free-boundary.

Each of these two types of molecules can regulate the progression of cells through the lineage by controlling their replication probabilities, through functional forms of and such as the following, In this case, is inhibited by with an EC50 value, or the half-maximal effective concentration, of and may achieve a maximal value of , while inhibits with an EC50 of and a thresholding value of . Hill exponents are represented by and .

Each of the unscaled cell lineage Eqs. (65–67) can be written in terms of , and generally aswhere and represents the corresponding right-hand side of each equation. Applying the transformation leads toBecause the growth of a tissue occurs at a significantly slower rate than that of a biochemical reaction, the morphogen system in Eqs. (69–74) can be evaluated at a quasi-steady state When the basal lamina resides along , the transformed boundary conditions in Eqs. (71–72) becomeThe Neumann boundary condition for a species along ,after the transformation takes the following form, which applies to both the morphogens and and also the cell densities .

To solve the convection equations describing the distributions of the cell densities in the tissue given by Eq. (78), a second-order upwind discretization was implemented in space along with a temporal discretization using a second-order TVD Runge–Kutta [23]. For the morphogen system in Eqs. (79–80), second-order central difference approximations are used in the - and -directions along with a multigrid solver for the reaction-diffusion equations. Along the dynamic boundary , a second-order upwind approximation is used in the -direction for terms with mixed partial derivatives. As is the case in solving for , 80 Gauss–Seidel iterations are used in the multigrid algorithm to solve for and . The system as a whole, with pressure and velocities solved using methods outlined in Sect. 2, maintains approximations in space and time that are at least second-order accurate.

First, we conduct numerical tests for the coupled system with the cell densities, pressure, dynamic boundary , and diffusive molecules and by assuming a rigid stem cell niche, a linear scaling of in Eq. (9), and nonuniform initial states of cell distributions and tissue morphology. Table 7 shows a spatiotemporal second-order accuracy for , , and in the numerical simulations for the stratified epithelia with a rigid niche. Because an exact solution cannot be constructed, errors are computed using the maximum-norm difference between successive approximations. Similar tests with a free-form stem cell niche also achieved a second-order accuracy for each of the system’s components in both time and space as evident in Table 8. For the tests with a free-form niche, an exponential form of is used to place more grid points in the -direction in a region near both near the interface and where sharp gradients of and form.

Table 7

Errors, orders of accuracy, and CPU times for calculation of epithelial growth and stratification with a rigid stem cell niche described in Eqs. (4–7, 65–72, 75–76) up to time

N	C_0		P		h		A		CPU (s)
	Error	Order	Error	Order	Error	Order	Error	Order
8	–	–	–	–	–	–	–	–	8.90e−2
16	1.66e−2	–	2.04e−5	–	1.66e−4	–	9.18e−4	–	1.38
32	6.22e−3	1.42	7.20e−6	1.50	3.55e−5	2.22	1.67e−4	2.45	24.2
64	1.54e−3	2.02	1.97e−6	1.88	7.53e−6	2.24	3.89e−5	2.10	396.0
128	3.84e−4	2.00	4.98e−7	1.97	1.67e−6	2.17	9.56e−6	2.03	7,499.5
256	9.51e−5	2.01	1.25e−7	1.99	3.90e−7	2.10	2.74e−6	1.80	142,867.4

A linear scaling is used for the transformation in Eqs. (8–10). is taken for time discretization. Initial conditions: , and . Chosen parameters are similar to those in [11]: , , , , , and . All hill exponents are chosen to be 2

Table 8

Errors, orders of accuracy, and CPU times for calculation of epithelial growth and stratification with a rigid stem cell niche described in Eqs. (4–7, 65–70, 73–76) up to time

N	C_0		P		h		A		CPU (s)
	Error	Order	Error	Order	Error	Order	Error	Order
8	–	–	–	–	–	–	–	–	9.50e−2
16	5.34e−2	–	7.17e−4	–	8.53e−5	–	1.22e−2	–	1.40
32	1.63e−2	1.71	1.87e−4	1.94	3.71e−5	1.20	3.24e−3	1.91	25.1
64	3.88e−3	2.08	4.68e−5	2.00	9.84e−6	1.91	8.19e−4	1.98	392.4
128	9.37e−4	2.05	1.16e−5	2.01	2.45e−6	2.01	2.07e−4	1.98	6,930.0
256	3.00e-4	1.64	2.89e−6	2.01	5.91e−7	2.05	5.35e−5	1.95	141,816.8

An exponential scaling is used for the transformation in Eqs. (8–10) with . is taken for time discretization. Initial conditions and parameters are given in Table 7

When beginning with an initial uniform distribution of stem cells and a distorted tissue morpohology, long-term behavior of the model with a rigid stem cell niche reaches an approximate steady state displaying stratification of cell layers in the -direction and homogeneity in the -direction. For the same initial conditions with a free-form stem cell niche, the tissue morpohology becomes further distorted and unstable fingers form in the tissue with stem cells accumulating at the tips. First, the tissue grows and increases in magnitude to time , and then the initial perturbation to amplifies over time as shown in Fig. 2. Plots of the distributions of each cell type in Figure 3 reveal that the stem cell niche and TA cell clusters localize at the tip of fingers where and concentrations are low. Thin tissue regions are, also, primarily dominated by TD cells, as low stem and TA cell quantities are present there.

Fig. 2

The dynamic interface at times for stratified epithelium simulations with a free-form stem cell niche. Initial conditions are , , and . A discretization size of is used along with . Parameters used are given in Table 7

Our method incorporating high-order central difference approximations and multigrid solvers enables fast and robust computation of stem cell tissue models in which biological processes occur at the molecular, cellular, and tissue level. This robust framework can also capture large distortion of tissue boundaries over long time scales by accurate and explicit tracking of the interfacial motion in the system.

An Expansion–Repression Model for the Imaginal Wing Disc in a Closed Geometry

The Drosophila imaginal wing disc serves as a popular model system for studying morphogen gradients and their role in long-range extracellular signaling in a developmental tissue. In particular, the morphogen decapentaplegic (Dpp), a bone morphogenetic protein (BMP) homolog, and its effects on wing vein placement and on growth of the wing disc are at the center of many recent experimental and modeling studies [3, 29, 34]. Dpp is produced along a thin region oriented along the dorsal-ventral axis and, as a result, its gradient forms along the anterior-posterior axis of the wing disc during the larval stages of development [55]. It has been shown experimentally that the gradient persists as the wing grows and maintains approximately the same length scale even as it assumes different sizes [59].

Recent work incorporates experiments and mathematical modeling to demonstrate that the morphogen Pentagone (Pent) diffuses throughout the disc and interacts extracellularly with Dpp in an expansion–repression mechanism to properly scale the Dpp gradient relative to the size of the disc [7, 26]. One of the models used to study this interaction, which does not explicitly include growth of the disc, is a two-component system consisting of a morphogen, , which serves as a repressor, such as Dpp, and an expander, , like Pent, in which the expander slows the degradation of the repressor while the repressor inhibits the production of the expander [7].

To account for the tissue growth, the equations for morphogens [7] need addition of convective terms, and the equations become Here, and are respective diffusion coefficients; and are respective maximal degradation rates; is the maximal production rate of ; and are corresponding half-maximal effective concentrations (or EC50s) for the feedback represented by a Hill function with as a Hill exponent. The function describes the production of Dpp in the disc with a smooth formin which is the maximal production rate.

To study the expansion–repression mechanism on a growing closed geometric domain, it is natural to use polar coordinates to represent the equations describing the imaginal wing disc. Similar to the morphogen systems in stratified epithelia, we assume that the time scale in which the morphogen reaches a steady state is significantly quicker than that of the tissue growth. Consequently, the morphogen system can be solved at a quasi-steady state, and, after applying the transformation, Eqs. (85–86) in polar coordinates become, No-flux boundary conditions along our dynamic boundary similar to the case without growth are imposed on the dynamic boundary, , with the following form, where or . To compute these transformed equations, we first use second-order central difference approximations in both the and directions. Since each equation is linear in one variable, we treat the resulting nonlinear system using a linear multigrid algorithm and a fixed point iterative solver.

To test the accuracy of our implementation, we first impose uniform proliferation by taking to be constant. Numerical simulations in which the disc is allowed to grow and then the nonlinear expansion–repression system is computed at the end of the simulation after growth has occurred, as opposed to at each time step, demonstrate a second-order accuracy in space and time as shown in Table 9. The observation in which order of three is obtained at larger in Table 9 suggest that the error may be mainly dominated by a dominance of error in the -direction, in which fourth-order approximations are used to compute the movement of and the internal pressure of the tissue.

Table 9

Errors, orders of accuracy, and CPU time for calculations of growth and the expansion–repression system on the imaginal wing disc given in Eqs. (35–37, 85–86) up to time up to time

N	P		h		[M]		CPU (s)
	Error	Order	Error	Order	Error	Order
16	–	–	–	–	–	–	0.571
32	1.20	–	4.73e−2	–	22.5	–	8.62
64	0.132	3.18	1.10e−2	2.11	22.5	0.82	140.4
128	1.70e−2	2.96	2.23e−3	2.30	3.02	2.07	2,595.9
256	3.52e−3	2.27	4.89e−4	2.19	0.374	3.02	48,841.2

A linear scaling is used for the transformation in Eqs. (38–40) along with . is taken for time discretization. Initial conditions: . Parameters used are adapted from the original model in [7]: , , , and . The Hill exponent is chosen to be 4

Experimental evidence suggest that the slope of the Dpp gradient prompts uniform proliferation of the wing disc as it grows [52]. Such regulation in growth in principle may be modeled via additional modifications in equations for . Because the focus of this study is numerical methods for tracking the moving boundaries, we instead impose a spatial distribution of influx of cells in such that the disc maintains its relative shape as it grows. In particular, we allow the tissue to grow until the anterior–posterior axis approximately doubles in length at . Figure 4a, b shows the distribution of on the disc at times and , respectively. To study how the length scales of the morphogen gradient compares, each of the scaled distributions are plotted on the unit disc using our linear transformation in Fig. 4c, d. It, then, can be seen that presence of an expander may increase the actual length scale of the gradient of as the disc grows in two dimensions, as similarly observed by one-dimensional models [7]. This notion is also evidenced by higher -values after growth has occurred. Whether in an open or closed geometry, our computational framework allows for the incorporation of diffusive morphogens and extracellular signaling, which enables tissue level models to account for crucial biological processes at the molecular level.

Fig. 4

Distributions of on the imaginal wing disc a prior to and b after growth of the disc. is also plotted on the unit square after scaled by c prior to and d after growth of the disc. A proliferation profile of is assumed. A discretization size of is used along with . Initial conditions and parameters used are given in Table 9 aside from , and

Conclusions

Spatial modeling of the stem cell niche and growth in developing and regenerative tissues consisting of multiple physical and biochemical processes is computationally challenging. In this paper, we have presented a continuum model of moving boundaries for tissue growth involving stem cells, cell lineages, and diffusive feedback molecules. In particular, we have developed a robust and accurate numerical method using a transformation technique to capture one or two dynamic boundaries in open geometries and growing closed geometries that are represented in polar coordinates. One important feature of the method is its capability of placing more grid points near dynamic boundaries that are tracked explicitly for increased resolution through a scaling by the transformation function. Another key element is an integration of a high-order finite difference approximation and a multigrid iterative solver for solving the transformed incompressibility equation for tissue.

This approach enables accurate computation of a moving boundary where surface tension exists due to cell-to-cell adhesion and where intertissue signaling molecules diffuse across the interface, leading to a second-order accuracy of the overall method in both space and time. Most importantly, such an approach enables the robust computation of undulated boundaries during tissue growth that are driven by the leakage of diffusive feedback molecules through explicitly and accurately tracked moving boundaries and the molecules’ regulation of stem cell fates. In particular, our method is robust in comparison to an iterative pseudospectral approach that is accurate but less robust in resolving deformed interfaces. Applications to developing stratified epithelia and the imaginal wing disc demonstrate that the overall numerical method is flexible in incorporating various biochemical processes (e.g. multiple morphogens) and capable of simulating long-term morphological tissue distortion.

In several instances, adaptive mesh refinement is used for the computation of interfacial motion computation in biological systems for necessary accurate boundary treatment [10, 14], but this approach comes with a difficulty in implementation and, often, the inclusion of pre-existing software packages [41, 62]. Our method allows for the accurate treatment of the boundary in a natural fashion by clustering grid points near dynamic interfaces using mathematical transformations to eliminate the necessity to include adaptive mesh refinement. However, with this approach comes the more difficult treatment of differential operators after transformation and also the constraint that the interface must be described as the graph of a function of or alone.

The modeling and computational framework presented can also be extended to other complex biological systems, such as the regenerative roles of stem cells in the liver and heart, or stem cell differentiation during tooth development. Importantly, each of the respective computational domains of these possible applications include interfaces that may be described as a graph of a function of one spatial variable and may contain more mechanistic effects, several lineages, and more complex downstream regulatory networks. The methods presented here can also be extended to geometries with non-simply connected domains, which may be relevant for tissue culture models in vitro in which intertissue communication by morphogens occurs, by considering separate pressure and interfacial variables for each simply connected subdomain. For efficient exploration of such complex models, it is important to refine the multigrid algorithms to account for variable coefficients in the tissue incompressibility equations such that the number of iterations during relaxations required for convergence becomes less sensitive to the size of the computational grid. More advanced temporal integrators, such as semi-implicit integration factor methods designed for stiff systems [45, 63], may be needed for better treatment of multiple time scales that are associated with many biochemical processes during tissue growth. Finally, it would be interesting to extend the model and computational framework to three spatial dimensions, providing a better understanding and additional insight into tissue morphologies and the roles of stem cells in growth in the third spatial direction that cannot be captured by two-dimensional models [27].