Evaluation of Diffusive Transport and Cellular Uptake of Nutrients in Tissue Engineered Constructs Using a Hybrid Discrete Mathematical Model

Tissue engineering systems for orthopedic tissues, such as articular cartilage, are often based on the use of biomaterial scaffolds that are seeded with cells and supplied with nutrients or growth factors. In such systems, relationships between the functional outcomes of the engineered tissue construct and aspects of the initial system design are not well known, suggesting the use of mathematical models as an additional tool for optimal system design. This study develops a reaction-diffusion model that quantitatively describes the competing effects of nutrient diffusion and the cellular uptake of nutrients in a closed bioreactor system consisting of a cell-seeded scaffold adjacent to a nutrient-rich bath. An off-lattice hybrid discrete modeling framework is employed in which the diffusion equation incorporates a loss term that accounts for absorption due to nutrient uptake by cells that are modeled individually. Numerical solutions are developed based on a discontinuous Galerkin finite element method with high order quadrature to accurately resolve fine-scale cellular effects. The resulting model is applied to demonstrate that the ability of cells to absorb nutrients over time is highly dependent on both the normal distance to the nutrient bath, as well as the nutrient uptake rate for individual cells.


Introduction
In tissue engineering applications for orthopedic tissues, such as articular cartilage or bone, a common bioreactor system consists of cells seeded in a biomaterial scaffold that is supplemented with nutrients or growth factors.In such systems, a diverse set of biophysical and physiological mechanisms interact to govern the formation of extracellular matrix that gives rise to the tissue's structure and resilient load bearing properties.The relationships between functional outcome measures for engineered tissue constructs and the design properties and conditions for the cell-scaffold system are not well understood.Trial-and-error approaches based on experiments involving a diverse set of cell-scaffold culture conditions can be time intensive and costly.The development of mathematical models for evolution of primary quantities of interest in bioreactor systems, as well as their incorporation into the design process, have the potential to accelerate the path to realization of optimal functional outcomes in engineered tissues; see [1] for a review of modeling approaches.
Interactions between the biological and physical mechanisms that govern tissue engineering outcomes occur at the microscopic scale of individual cells.On this scale, cell-cell and cell-matrix interactions are typically time varying, nonlinear functions of the variables involved.While, in theory, it may be possible to systematically quantify ensemble effects based on multiscale modeling techniques, such as numerical homogenization [2,3], this approach typically requires a priori specification of the constitutive forms of relations among dependent variables in the system.These constitutive relations may not be easily identified and, overall, such techniques have seen little application in modeling orthopedic tissue engineering in cell-scaffold systems.Several models for cartilage tissue engineering are based on representing the evolving system as a coupled set of biological reactions idealized as a system of ordinary differential equations [4,5].Other approaches are spatio-temporal and based on continuum representations via partial differential equations (PDEs) within the framework of reaction-diffusion models [6][7][8][9][10] or mixture theories [11].It should be noted that these studies involve the representation and solution of a forward model, and ultimately, the design of tissue engineered constructs necessitates the integration of such models into with optimality criteria; see e.g., [12].
In modeling other types of tissues, an alternate approach involves the explicit representation of individual cells, with the total number of cells being sufficient to simulate and identify tissue-scale phenomena of interest in the system [13].In lattice-based models, each lattice cell on a regular grid represents a region containing one constituent (e.g., cell, scaffold, extracellular matrix) in the system, and the total energy is formulated on the grid and evolved via an iterative algorithm.In contrast, off-lattice models do not employ a grid, but represent each cell as a subdomain within the system of macroscopic governing equations that model the evolution of the engineered tissue construct.In this study, we employ the off-lattice approach based on a hybrid discrete modeling framework [14][15][16].Hybrid discrete models [17][18][19][20][21] enable coupling of PDE continuum models for variables, such as nutrient concentration, with models that treat cells as individual discrete entities for which biological or physiological dynamics can be separately formulated.
We focus on two primary mechanisms in the cell-scaffold system, which are the diffusive transport of nutrients and the loss of nutrients due to cellular uptake.The latter mechanism is assumed to be required for the maintenance of cell viability and to sustain biosynthesis and cross-linking of extracellular matrix proteins.In the associated hybrid discrete reaction-diffusion model, each individual cell has a small region of influence in which nutrients are absorbed from the background medium.Robust numerical solutions of the associated governing equations are employed and based on the use of discontinuous Galerkin finite elements in space and a Crank-Nicolson finite difference method in time.Simulations are performed to evaluate the competition for nutrients between cells located near and far from a nutrient source, noting that tissue regeneration is known to occur earlier and faster in peripheral regions adjacent to a nutrient-rich bath [22].

Model Description
We use a hybrid discrete modeling framework [14][15][16] in which a continuum model for nutrient diffusion with a concentration of u(x, t) is augmented with source terms accounting for individual cells that remove nutrients as the system evolves in time.The following nonlinear PDE model is considered in a two-dimensional domain, Ω, and for t > 0: Here, f denotes the loss term that models the uptake of nutrients from the domain, Ω, due to the presence of N c individual cells, each centered at 1) is of the hybrid-discrete type, treating each cell as a discrete entity that absorbs nutrients according to a prescribed nutrient dependent function, g(u), where u b in Equation (1) represents a background nutrient concentration.For each individual cell, this absorption occurs over a region with spatial extent specified by a support function, δ , with an associated characteristic length scale, .We associate the length scale, , with a specified cell radius that, effectively, determines the volume fraction of cells in the engineered tissue.The positive constants b i (i = 1, . . ., N c ) control the nutrient absorption rate and may vary for each cell.
Our model assumes a two-dimensional rectangular domain , where 1 mm is taken as the unit of length.The domain, Ω, is sub-divided into two subdomains consisting of one region on the left representing a nutrient-rich fluid bath 1).The support function, δ , for individual cells is assumed to cover a two-dimensional disk, C i , (i = 1, . . ., N c ), of radius = 5 µm, and centered at the cell centers, x c i , (i = 1, . . ., N c ).In constructing the model, cell centers are located in Ω 2 in a random, non-overlapping, manner.Once all cells have been placed in the domain, we define the cell volume fraction, Φ c , as the ratio of the total area occupied by cells divided by the area of Ω 2 .Each cell is also assumed to remain stationary within Ω 2 and to maintain its initial size throughout time.We assume an initial nutrient concentration profile that is a positive constant on the left half of the domain, Ω 1 , and zero on the right half of the domain, Ω 2 (Figure 2).No flux boundary conditions ∂u/∂n = 0 are prescribed on the entire boundary, ∂Ω, corresponding to an assumption that there is a fixed amount of total nutrients available in the bioreactor.Lastly, in modeling diffusive nutrient transport, we assume a globally constant diffusivity, D i , in each subdomain, Ω i (i = 1, 2), and that, in region Ω 2 , the diffusivity D 2 is the same in both the cellular and extracellular regions.In the left half, Ω 1 , the normalized nutrient concentration is assumed to be u = 0.2, and on the right half, Ω 2 , the initial nutrient concentration is taken to be zero.
Detailed modeling of the complex intracellular dynamics that may be involved in nutrient uptake by individual cells is beyond the scope of this study.Instead, to account for the effects of nutrient uptake by individual cells, we consider a simplified representation of the cell supported function, δ (x − x c i ), in Equation (1) that represents the spatial component of the cell model.Specifically, for i = 1, . . ., N c , we assume a uniform cell model: (Figure 3) This model assumes that the absorption of nutrients in the bioreactor system occurs in a spatially uniform manner within all regions, C i , that represent the individual cells in Ω 2 .Lastly, a model is needed for the function, g(u), in Equation ( 1) that governs the cell-mediated regulation of nutrient uptake within each individual cell.Our model is formulated by multiplying the cell supported function, δ , by a normalized nonlinear nutrient dependent loss term, g(u) [23], that is assumed to have the form: and g(u) = 0 for u < u * (Figure 4).Here, u * is prescribed a priori and represents a minimum nutrient concentration for cell viability, below which the cell ceases to absorb nutrients.The parameter λ > 0 regulates the rate of nutrient absorption as a function of local nutrient concentration within an exponential model (3).We non-dimensionalize our governing Equation (1) in a standard way by dividing u(x, t) by the background concentration, u b , taking a = 1 mm as the unit of length and taking a 2 D as the unit of time, where D = max {k=1,2} D k .Transforming to non-dimensional variables, we have: Based on this transformation, the following non-dimensional version of Equation ( 1) is obtained: where bi = b i a 2 D , ¯ = a and Dk = D k D (k = 1, 2).In the subsequent sections of this paper, for simplicity of notation, we drop the bar symbols.

Numerical Methods
Since our governing equations are nonlinear, inhomogeneous at the fine scale and do not possess analytical solutions, we require numerical methods that are both stable and accurate.Consequently, we employ high order discontinuous Galerkin (DG) finite element numerical schemes.Since finite elements generate a variational formulation of Equation ( 5), we can readily handle the effects of the cell supported functions, δ , via numerical quadrature and adequate levels of mesh refinement (see Figure 5).High order solution approximations can be obtained by using higher degree basis functions.Specifically, due to its local nature, the DG finite element method allows us to express and work with our quantities on single mesh elements without the requirement of additional globally prescribed conditions.Another attractive feature of the DG method is the ease with which it accounts for various types of boundary conditions by incorporating them into the weak formulation.In contrast, regular finite elements typically require restrictions on the finite element space, e.g., to enforce essential boundary conditions.Lastly, in the context of adaptive mesh refinement, the DG method produces meshes with fewer degrees of freedom, as it allows hanging nodes, whereas continuous finite element methods may require unnecessary refinement in order to preserve mesh conformity.We now briefly describe the specific numerical method utilized to solve the governing equations in this study.

Figure 5.
Numerical treatment of cellular effects in the hybrid discrete model.(A) High order quadrature with 37 points and based on a thirteenth-order Gaussian quadrature rule, is used to adequately resolve the effects of the cell supported δ functions in the element integrals; (B) Illustration of a typical uniform mesh taken from a portion of the subdomain, Ω

B
There is an assortment of DG variants that can be found in the literature.For a survey on DG methods for elliptic problems, see [24][25][26][27], and for parabolic problems, see [28,29].In this study, we have used the symmetric interior penalty discontinuous Galerkin finite element (SIP-DGFE) method for the spatial discretization and the Crank-Nicholson finite difference scheme for the time discretization of (1).The SIP-DGFE method produces, in primitive variable form, symmetric operators that can provide greater flexibility during implementation.The following fully discrete numerical scheme corresponding to (1) computes the discrete approximation, U m , of the nutrient concentration function u at t = t m : Let 0 = t 0 < t 1 < . . ., t M = T be a partition of [0, T ], and Here, , where V h is the DG finite dimensional space [25] consisting of polynomial element basis functions, τ m is the time step at the m-th time iteration and α h (•, •) is the DG bilinear form corresponding to the Laplacian operator in variational form with prescribed boundary conditions (BC) [26].It should be noted that, in our model, the fact that we have a priori knowledge of the boundaries of the subdomains is crucial in choosing an appropriate triangulation for the entire domain, such that the edge of the triangles follows the internal boundary that separates regions Ω 1 and Ω 2 .Specifically, within the context of DG finite elements, the bilinear form can naturally incorporate the known discontinuity in the diffusion coefficient for such a triangulation.The q-th-order DG basis polynomials, e.g., q ∈ {1, 2, 3, 4}, provide accuracy that is O(h q+1 ) for the L 2 norm of the spatial error at each time step, with h being the mesh discretization parameter.In this study, we have set q = 3, which, in theory, gives us an L 2 -spatial error that is O(h 4 ).Time discretization using the Crank-Nicolson finite difference scheme provides second order accuracy in time, i.e., O(τ 2 ), where τ = max i τ i .Hence, overall, we expect to have a global L 2 -error of O(h q+1 + τ 2 ).We note here that standard convergence tests have been performed to gauge the aforementioned order and found that our schemes performed as expected.Having the formulation in primitive variable form makes it ideal to exploit the symmetry and positive definiteness of the associated operator.We have used an implicit-explicit inner iteration scheme in order to treat the nonlinearities present in Equation (6).A conjugate gradient solver with multigrid preconditioning (PCG) is used to solve the resulting linearized algebraic system arising from our numerical discretization.The code is also capable of performing mesh refinement and coarsening, as needed, by employing a marking strategy based on inverse estimates, described in [30][31][32].

Results
The model and numerical techniques described above were used to evaluate interactions between diffusive transport of nutrients and nutrient loss due to cellular uptake in the system represented in Figure 1.Our focus was to quantitatively evaluate within region Ω 2 the effect of distance from the interface to region Ω 1 on the ability of cells to absorb nutrients over time.An averaging operator is introduced to quantify the amount of nutrients present in each cell at a given time, t.Specifically, for any function, κ ∈ L 1 (C i ), define the operator: which computes an average value of function κ over each disc, C i , i = 1, . . ., N c .Based on Equation (7), at a specific time, t m , we calculate the nutrient being consumed by a cell, i, as: Hence, the total amount of nutrients consumed by cell i during the entire simulation is given by: where M is the total number of time steps in the finite element numerical scheme.In all simulations performed, the non-dimensional values of the diffusivities were set as D 1 = 1 (in Ω 1 ) and D 2 = 0.5 (in Ω 2 ), and the total time was chosen as T = 3.5.In non-dimensional units, parameters for the cell model were taken as = 0.005 and λ = 10 (see Equation ( 3)).The non-dimensional value of the cell absorption rate was then varied among the non-dimensional values b i = b = 0.5, 1, 2.5, 5. Overall, the total amount of nutrient absorbed per cell, Λ i , varied with both spatial location and as the nutrient absorption rate, b, was changed.The total amount of nutrient consumed by each cell is shown as a function of spatial location within region Ω 2 in Figure 6 for the case b = 0.5.We observe that there is variation among the total absorption values in both spatial directions (X and Y ), but that cells closer to the interface (X near zero) with the nutrient bath region, Ω 1 , consume substantially more nutrients.For this particular case, cells furthest away from the interface (X near one) absorb almost ten times fewer nutrients than the cells that are closest to the interface with the nutrient source.We also examined the effects on the total nutrient consumed by each cell as values of the nutrient absorption rate, b, were varied.The cases for b = 0.5, 1, 2.5, 5 were simulated and are plotted in Figure 7 as a function of the normal distance (X) to the interface between Ω 1 and Ω 2 .These results clearly illustrate a strong dependence of the total amount of nutrients consumed per cell on both the normal distance to the nutrient source (X direction) and the nutrient absorption rate, b.For each fixed value of b, data points for total nutrient consumed by all cells were curve fit to a fourth degree polynomial obtained using MATLAB's "cftool" package (see Figure 7).A comparison of the resulting polynomial fits indicates that cells distant from the source of nutrients may have limited access to nutrients, independent of cellular absorption rate, due to their inability to compete with cells closer to the supply for access to nutrients.

Discussion and Conclusions
In this study, a reaction-diffusion model was developed to study the competing effects of diffusive nutrient transport and cellular uptake of nutrients in a closed bioreactor system comprised of a cell-seeded scaffold adjacent to a nutrient-rich bath.The model was formulated within the framework of off-lattice hybrid discrete models [14][15][16] that couple cellular effects to PDE models for nutrient diffusion via a loss term that accounts for nutrient absorption by each individual cell.Numerical solutions of the resulting governing equations were developed based on the use of a discontinuous Galerkin finite element method with high order quadrature being used to accurately resolve fine-scale variations due to the microscopic effects of individual cells.The resulting numerical models were applied to demonstrate that the normal distance of cells from the source of nutrients, as well as the nutrient uptake rate for individual cells both strongly influence the ability of cells to absorb nutrients over time.These findings are consistent with prior experimental studies demonstrating that regions of engineered cartilage constructed nearer to nutrient access points exhibit tissue growth faster than other regions of the tissue [8,22].
The model of an engineered tissue construct presented in this study is based on several simplifying assumptions.These include the restriction of the model to two spatial dimensions and the assumption that intracellular and extracellular diffusivities were constant and equal.Furthermore, only two primary mechanisms, nutrient transport due to diffusion and nutrient loss due to cellular uptake, were considered.As such, the focus of this study was to demonstrate the potential benefits of a hybrid discrete modeling approach in quantifying the collective effects, at the macroscopic scale, of cells whose dynamics are specified via a distinct set of models formulated at a finer scale.A more detailed and accurate model of cell-scaffold systems for orthopedic tissue engineering would require the incorporation of many additional mechanisms into the models, including the balance of mass for bound and unbound system constituents, the balance of momentum, intracellular dynamics and cell proliferation.In addition, recent experimental studies for cartilage tissue engineered in hydrogels demonstrated that initial pore size can strongly affect functional tissue outcomes [33], suggesting that mixture modeling approaches [4,11] may be an appropriate framework for the development of more detailed models.While the modeling framework proposed in the current study extends to model formulations that are more detailed and in three dimensions, substantial increases in the computational cost of the associated numerical solutions are to be expected and the amenability of both the modeling approaches and numerical methods to parallelization and mesh adaptivity are important factors to consider in future model extensions.

Figure 2 .
Figure2.Illustration of the initial nutrient concentration profile on the domain, Ω.In the left half, Ω 1 , the normalized nutrient concentration is assumed to be u = 0.2, and on the right half, Ω 2 , the initial nutrient concentration is taken to be zero.

Figure 4 .
Figure 4. Illustration of the nutrient-dependent model for nutrient loss g(u) (3) as a function of non-dimensional nutrient concentration in the cell model.The case with λ = 10/u 2 b and minimum nutrient concentration u * /u b = 0.05 is shown, where u b is a background nutrient concentration used in the non-dimensionalization.

Figure 6 .
Figure 6.Illustration of the total nutrient consumed by each cell in the cell-seeded scaffold (Ω 2 ) for the case b = 0.5.

Figure 7 .
Figure 7.The effects of normal distance (X direction) to the source of nutrients and nutrient absorption rate b on the total amount of nutrient consumed by each cell.Fourth degree polynomials were fitted to the values of total nutrient consumed by each cell in the cases b = 5 (red), b = 2.5 (black), b = 1 (green) and b = 0.5 (purple).