Analytic Models of Oxygen and Nutrient Diffusion,Metabolism Dynamics,and Architecture Optimization in Three-Dimensional Tissue Constructs with Applications and Insights in Cerebral Organoids

Quasi-steady-state solution in linear one-dimensional cases

Three-dimensional shapes have often been modeled with one-dimensional (1D) solutions in applications such as thermodynamics and diffusion, including in common configurations of planar, cylindrical, and spherical constructs,^38–42 and numerically solved 1D diffusion models have been shown to be an accurate and valid simplification for many types of 3D tissue constructs.^43,44 One-dimensional analytic solutions are also applicable to diffusion in tissue constructs, particularly in constructs where the majority of diffusion occurs in a uniform direction (designated as distance x, see Fig. 1) or where diffusion occurs through a primary surface such as the top layer of a planar or slab-like construct. One-dimensional solutions are also valid when diffusion occurs through the outer surface of a cylinder along a radial axis in cylindrical coordinates or through the outer surface of a spherical construct in spherical coordinates along the radial axis (designated as radial distance r, see Fig. 1). Solutions with circular or spherical symmetry are producible by altering the coordinate frame and adjusting boundary conditions appropriately for the diffusion system, thereby accurately depicting diffusion and metabolism for specific forms of 1D, 2D, or 3D constructs.

OPEN IN VIEWER

FIG. 1.

Views of three-dimensional (3D) tissue constructs in a culture well. The thickness of the planar construct is T and the diameter of the sphere is 2R.

Diffusion can be modeled using physical laws originally described by Adolf Fick, depicting a changing flux over distance and conservation of mass over distance and time. Fick had a passion for mathematics and physics, and although he chose to become a physician, he continually sought to apply the concepts of mathematics to the study of the human body. Fick's model is based on similar forms of equations for transfer of heat and energy, as described by Joseph Fourier and Isaac Newton, respectively, being derived from the conservation of mass within a defined region, wherein the change in concentration (C) per time (t) is equal to the change in flux (J) over distance (x), written as \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$ \frac { \partial C } { \partial t } + \frac { \partial J } { \partial x } = 0$$ \end{document} . By substituting Fick's first law for flux \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$J \,=\, - D \frac { \partial C } { \partial x } $$ \end{document} into the conservation of mass equation (where D is diffusivity), the classic form of the diffusion equation is obtained: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \frac { \partial C } { \partial t } = D \frac { \partial^2C } { \partial x^2 } \tag { 1 } \end{align*} \end{document}

To describe exact solutions to this equation, the 1D case of distance (x) into a tissue construct is first considered here (Fig. 1). The solutions to Eq. (1) can be described in steady state and nonsteady state and can also be given over the infinite, semi-infinite, or finite spatial domains. By maximizing the tissue construct thickness such that complete consumption of nutrient occurs at the deepest point in the construct, certain difficulties in the mathematics can be surmounted. A useful quasi-steady state solution (where the change in concentration of a nutrient in time is constant) can be found for a tissue construct that has characteristics of both diffusion and metabolism by first adding the change in concentration over time to a metabolic component called φ, representing the consumption rate of the nutrient by the cells in culture (in units of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$ \frac { mol } { L \cdot s } $$ \end{document} herein). By then letting the diffusional component reach quasi-steady state \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\left( \frac { \partial C } { \partial t } + \varphi = \varphi \right)$$ \end{document} , the equation becomes \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}\varphi = D \frac { \partial^2C } { \partial x^2 } \tag { 2 } \end{align*} \end{document}

It is reasonable to treat φ as a constant (zero order kinetics) when cells are no longer dividing and are in a stable metabolic state, although in many conditions there are often phases of growth, differentiation, and consumption that may fluctuate over time. By defining C_o as the concentration of the nutrient in the media at the surface of the hydrogel, the surface of the hydrogel as x = 0, and T as the length through the hydrogel from the surface to the floor of the construct (Fig. 1), the maximal thickness of tissue construct can be solved by integrating twice. Certain boundary conditions are applied in this process, depending on the construct design. With an impermeable floor at the base of the construct (x = T), the change in concentration ceases, given as \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$ \frac { \partial C ( x = T , t ) } { \partial x } = 0$$ \end{document} . Furthermore, with the assumption that the concentration at the surface of the construct remains constant, given by \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$C ( x = 0 , t ) = C_o$$ \end{document} , the integration of Eq. (2) becomes \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}C ( x ) = \frac { \varphi x^2 } { 2D } - \frac { \varphi Tx } { D } + C_o \tag { 3 } \end{align*} \end{document}

Substituting C(x) = 0 at x = T represents the distance at which the nutrient is fully consumed, which then further defines the maximum thickness of diffusion in a cellularized hydrogel: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}T_ { max } = \sqrt \frac { 2C_oD } { \varphi } \tag { 4 } \end{align*} \end{document}

If diffusion is allowed from both below and above the construct (such as with planar constructs in a suspended well in the media rather than seated on the floor ³⁴ ), the T_max is simply doubled.

Quasi-steady-state solution in radial cases

Radially symmetric constructs may also be modeled with 1D solutions by substituting spherical or cylindrical coordinates into the diffusion equation. In the spherical case, just as in aforementioned conditions, it is assumed that density and diffusion coefficients will be constant throughout the construct and the dimensions of the sphere are constant without swelling or shrinking. By substituting spherical coordinates into the 3D diffusion equation \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$ \frac { \partial C } { \partial t } = D \Delta C$$ \end{document} (where Δ is the LaPlacian operator, or the divergence of the gradient of C) and by implementing azimuthal and polar symmetry around the sphere (meaning the change in concentration with respect to θ or φ is zero, as would be expected in an ideal sphere), diffusion in a sphere of radius R along the radial axis r can be written as follows^38,42,45: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \frac { \partial C } { \partial t } = \frac { 1 } { r^2 } \frac { \partial } { \partial r } \left( r^2 D \frac { \partial C } { \partial r } \right) \tag { 5 } \end{align*} \end{document}

In fact, a simpler way of writing this equation for any dimensionality s is ⁴⁶ \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \frac { \partial C } { \partial t } = \frac { 1 } { r^ { s - 1 } } \frac { \partial } { \partial r } \left( r^ { s - 1 } D \frac { \partial C } { \partial r } \right) \tag { 6 } \end{align*} \end{document}

where s = 1 describes the 1D axis, s = 2 describes a cylinder where diffusion at the end surfaces is negligible and occurs primarily along a circular radius, and s = 3 describes a sphere. Assuming a constant nutrient consumption rate and isotropic diffusivity sets up the spherical equation similar to Eq. (2): \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}\varphi = \frac { D } { r^2 } \frac { \partial } { \partial r } \left( r^2 \frac { \partial C } { \partial r } \right) \tag { 7 } \end{align*} \end{document}

Integrating twice produces \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}C ( r ) = \frac { \varphi r^2 } { 6D } - \frac { c_1 } { r } + c_2 \tag { 8 } \end{align*} \end{document}

Boundary conditions are produced by the condition of spherical symmetry, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$ \frac { \partial C ( r = 0 ) } { \partial r } = 0$$ \end{document} , giving c₁ = 0, and from the initial concentration at the outer edge of the sphere (at point R), the concentration C(r = R,t) = C_o. Then c₂ can be solved as \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$c_2 = C_o - \frac { \varphi R^2 } { 6D } $$ \end{document} , giving a general solution for concentration along r in the construct: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}C ( r ) = \frac { \varphi } { 6D } ( r^2\, -\, R^2 ) + C_o \tag { 9 } \end{align*} \end{document}

The maximal depth occurs when C = 0 at r = 0, which can then be solved to find maximal R for a given φ: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}R_ { max } = \sqrt \frac { 6C_oD } { \varphi } \tag { 10 } \end{align*} \end{document}

The R_max thus represents the maximal diffusion distance from a cell to surface source of nutrient, and this has been used, for example, as the maximal viable radius in numerical models of spherical tumor growth.^47–49 By the same methods, the general cylindrical solution, where diffusion occurs through the surface curvature and not at the ends, is of the form \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}C ( r ) = \frac { \varphi r^2 } { 4D } - c_1 \ln r + c_2 \tag { 11 } \end{align*} \end{document}

Under the described boundary conditions, the general solution and maximum radius of a cylindrical construct then become \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}C ( r ) = \frac { \varphi } { 4D } ( r^2 - R^2 ) + C_o \tag { 12 } \end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}R_ { max } = \sqrt { \frac { 4C_oD } { \varphi } } \tag { 13 } \end{align*} \end{document}

A comparison of multidimensional steady-state solutions with the same parameters is shown in Figure 2.

OPEN IN VIEWER

FIG. 2.

Comparison of metabolic steady-state profiles of oxygen for thickness (T) or radius (R) in mm in 1D, 2D, and 3D models, each with identical values for φ, D, and C_o (with an inverted spatial axis for the radial cases for comparison).

General properties of a tissue diffusion system

The value of the diffusivity coefficient (D) will depend on the characteristics of both the diffusant and the medium through which the molecule diffuses. These characteristics may include nutrient properties like size, charge, and mass ⁵⁰ ; likewise, relevant medium properties may include variations in material density, cellular density, and temperature, all of which can be modeled in a nonlinear manner. Nevertheless, it may generally be assumed that the diffusivity coefficient remains a constant parameter in a tissue construct of uniform material density, meaning it is independent of factors such as changing concentration of solutes or cellular demand.^43,51 The models presented in this study also do not include swelling or shrinking of cells, tissues, or hydrogel material or other poroelastic effects, and the assumption of homogenous and isotropic tissue allows derivation of solutions that apply broadly to engineering tissue constructs. For various small molecules in hydrogel polymers, the diffusion coefficient is typically around 10⁻⁹ to 10⁻¹⁰ m²/s, with small molecules like oxygen typically closer to 10⁻⁹ and slightly larger molecules like glucose closer to 10⁻¹⁰ m²/s (different units of measure in prior studies are normalized here).^44,51–60 In hydrogels, diffusivity has been reported to be about 50% to 85% of that found in solution, and diffusivity in cartilage has been reported to be about 40% of that found in solution, with diffusion coefficients for small molecules in a similar range of 2 × 10⁻¹⁰ to 8 × 10⁻¹⁰ m²/s.^43,56,61,62 Larger proteins like albumin or glycosaminoglycans generally have a smaller diffusion coefficient in the range of 10⁻¹⁰ to 10⁻¹¹ m²/s.^62–64

The initial concentrations (C_o) of many types of nutrients are generally provided in the media solution, and culture media may vary considerably in its nutrient content, containing physiologic concentrations of glucose or several times higher concentrations to compensate for intermittent media changes and lack of continual perfusion, or to study certain conditions and cell types (see Table 1). In cell cultures, standard values of diffused oxygen and nutrient concentrations in gas and liquid phases surrounding the tissue construct may generally be assumed, especially when culture media is well mixed, or also, in the case of gasses diffusing from outside the media, when the layer of media over the tissue construct is thin or equilibrated for long times. The typical value of oxygen concentration in various depths of solution under standard atmospheric pressure and at 37°C has been reported to be between 0.20 to 0.22 mM.^{44,57,65–68} For greater depths of stagnant media over a construct, the concentration of oxygen at the surface of the tissue construct (which here, as with other nutrients at a construct surface, is called C_o) can be approximated using Fick's first law of diffusion through overlying liquid media, which, as before, states that flux is equal to the diffusion coefficient multiplied by the change in concentration over the change in distance. Using a thin film assumption, this can be written as \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}J = \frac { - D } { l } \left( \frac { P } { H } - C_o \right) \tag { 14 { \rm a } } \end{align*} \end{document}

where J represents the net flux of oxygen in the x-direction, l is the length or thickness of the media over the construct surface (e.g., in length units that match those used in the diffusion coefficient), and the oxygen concentration at the media surface is given by P/H where P is the partial pressure of oxygen (e.g., in units of atm) and H is Henry's constant for the media solution (e.g., in units of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$ \frac { L \cdot atm } { mol } $$ \end{document} ). Assuming quasi-steady-state conditions, where φ is the metabolic rate of the cellular construct (in units of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$ \frac { mol } { L \cdot s } $$ \end{document} ), substituting the flux from Eq. 14a into Eq. 2, and solving for concentration through the liquid media then gives \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}C = \frac { P } { H } - \frac { \varphi lx } { D } \tag { 14 { \rm b } } \end{align*} \end{document}

The value of C_o is found when x = l, and D is the diffusion coefficient of the gas in the surrounding liquid media, typically around 3 × 10^–9 m²/s for oxygen at physiologic temperature and pressure.^43,44

The rate of oxygen consumption ( \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$m_{o_2}$$ \end{document} ) of neuronal cells in brain across species has been reported to be ∼1 × 10⁻¹⁷ to \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$1 \times 10^ { - 15 } \frac { mol } { cell \cdot s } $$ \end{document} ^69–71 (Table 2). The rate of glucose consumption (m_gluc) in neuronal cells has been reported to average from about 8 × 10^–17 to \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$4 \times 10^ { - 16 } \frac { mol } { cell \cdot s } $$ \end{document} across a variety of rodent and primate species, with a lower range of 1 × 10^–17 to \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$2 \times 10^ { -17 } \frac { mol } { cell \cdot s } $$ \end{document} for cerebellar neurons and a higher range of 2 × 10^–16 to \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$4 \times 10^ { - 16 } \frac { mol } { cell \cdot s } $$ \end{document} for cortical neurons (Table 3). In humans, cerebellar neurons consume glucose at an average rate of about \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$1 \times 10^ { - 17 } \frac { mol } { cell \cdot s } $$ \end{document} while cortical neurons consume glucose at an average rate of about \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$2 \times 10^ { - 16 } \frac { mol } { cell \cdot s } $$ \end{document} . ⁶⁹ The metabolic rate of whole brain among species varies as a function of brain size and neuronal density, but the metabolic rate per cell is remarkably similar among species, as has been noted previously. ⁶⁹ These rates can be adapted for the density of cells within a tissue construct by using the formula φ = mρ, where m is the metabolic rate for the nutrient per each cell (e.g., in units of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$ \frac { mol } { cell \cdot s } $$ \end{document} ) and ρ is the average density of cells in the hydrogel or tissue (e.g., in units of cells/L). This formula is based on the premise that cells are homogeneously distributed or metabolic consumption values are similar throughout the construct, a premise which is generally valid for relatively uniform structures like neurospheres and embryoid bodies. ⁷² Thus, as an example, for a construct with 20,000 neuronal cells in 25 μL of hydrogel, ρ = 8.0 × 10⁸ cells/L, and the final value of φ for the construct becomes about 9.8 × 10^–9 to \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$8.8 \times 10^ { -7 } \frac { mol \ O_2 } { L \cdot s } $$ \end{document} for oxygen. Of note, the metabolic consumption by the entire construct can also be found by multiplying φ by the volume of the 3D tissue construct (e.g., \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$V_c \,=\, \frac { 4 } { 3 } \pi r^3$$ \end{document} for a sphere), which can be useful in determining the perfusion needed or the amount and frequency of media changes for nutrients supplied by the media.

The minimal and maximal thickness of hydrogel can then be estimated by these parameters and formulas. If minimizing metabolic parameter values are used—that is, if neuronal cells described above consume oxygen at an upper range of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$8.8 \times 10^ { - 7 } \frac { mol \ O_2 } { L \cdot s } $$ \end{document} and the diffusion coefficient is at the lower range of 10⁻¹⁰ m²/s—then the viable thickness of tissue is predicted to be only about 250 μm. Conversely, if the maximizing parameters are used, meaning oxygen consumption of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$9.8 \times 10^ { - 9 } \frac { mol \ O_2 } { L \cdot s } $$ \end{document} and diffusion coefficient of 10^–9 m²/s, then the maximum thickness of viable tissue under these idealized culture conditions could be as much as 16 mm assuming steady state with unconstrained flux at the construct surface. For the average human neuron metabolic rate of oxygen and at an intermediate diffusivity, the viable slab thickness is predicted to be 1.4 mm, with a range of 0.4–2 mm for minimal and maximal diffusivities, respectively, or thicker for more quiescent cells, meaning that the viability of neurons in 3D tissue depends significantly on the neural activity and the proximity and diffusivity of nutrients in tissue. This thickness limit for cell viability has been observed experimentally in layered hydrogel constructs, with differentiated human neurons with single surface layer diffusion, which were capable of being up to 2-mm thick. ³² Similarly, analysis of the spherical equation for a construct with 20,000 neurons in 25 μL of hydrogel provides corroborating results, with the viable diameter being 0.9–4.1 mm for minimal and maximal diffusivities, respectively. This means that a neuronal sphere with normal spontaneous neural activity and moderate metabolic activity and a diameter of about 4 mm should just be viable at this cell density, and it is interesting to note that this was the size limit reported in cerebral organoids made from pluripotent stem cells differentiated within spherical hydrogels, although the cellular density in these organoids was unknown. ³¹

Diffusion models (steady and nonsteady states) in the infinite domain

The diffusion equation can be solved through several well-known methods if it is assumed that diffusion occurs through an infinite spatial distance. Because these solutions are only applicable over an infinite spatial range, however, they are not particularly helpful for modeling diffusion and metabolism in finite constructs, but they are helpful for understanding the general concepts of diffusion solutions. Solution methods include the use of LaPlace transforms or, through the strategy of Boltzmann's transformation, substituting \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$ \mu {{def} \atop} \frac {x} {\sqrt {4Dt}} $$ \end{document} into Eq. 1 to make the partial differential equation (PDE) into an ordinary differential equation (ODE) ⁴² and solving through a series of steps with desired initial conditions. \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}- 2 \mu \frac { \partial C } { \partial \mu } = \frac { \partial^2C } { \partial \mu^2 } \tag { 15 } \end{align*} \end{document}

When integrated twice with the initial condition C(x,t = 0) = C_o over the semi-infinite distance ( \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$0 < x < \infty$$ \end{document} , the result is^42.52,73 \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}C ( \mu ) \,=\, C_o [ 1 - { \rm erf} ( \mu ) ] \tag{16}\end{align*} \end{document}

where \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} { \rm erf } ( \mu ) = \frac { 2 } { \sqrt { \pi } } \int^ { \mu } _0 e^ { - \vartheta^ { 2 } } d \vartheta \tag { 17 } \end{align*} \end{document}

wherein, the variable ϑ acts as a place holder for μ over the integration. It can be verified that this is indeed a solution by substituting back into the original Eq. 1. Finally, substituting back for μ gives the semi-infinite solution in terms of x \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}C ( x , t ) = C_o \left[ 1 - { \rm erf } \left( \frac { x } { \sqrt { 4Dt } } \right) \right] \tag { 18 } \end{align*} \end{document}

Although the error function (erf) appears complex, it can be solved by many types of programs, including simply typing “erf( )” into the website google.com. This provides a solution in the infinite spatial domain for diffusants that are not consumed by the cells and held constant in the surrounding solution [the condition (0, t) = C_o]. It can be seen here that the value \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\mu = \frac { x } { \sqrt { 4Dt } } $$ \end{document} is a dimensionless parameter, which allows comparison of diffusion behavior across different times and diffusivities. The time required for the concentration to reach a given value at a given point is proportional to the square of the distance into the construct and inversely proportional to the diffusivity. Values for desired time (t) and diffusivity (D) can be inserted and the concentration profile solved and graphed (Fig. 3e, f). Under the above assumption that the amount of oxygen and nutrient outside the hydrogel remains constant, this model shows that the concentration of diffusant reaches 50% of C_o through 1 mm of hydrogel in about 20 min and 90% of C_o in ∼8.8 h when D = 10⁻⁹ m²/s. This is in accordance with other reports of an 8-h diffusion time for small molecules over a distance of 1 mm in hydrogel. ⁷⁴

OPEN IN VIEWER

FIG. 3.

Concentration profiles in a semi-infinite domain of constant diffusivity as a function of distance x (depth through hydrogel in mm) and time t (in seconds). Having assumed a finite amount of diffusant (like glucose in a culture well, where D = 10⁻¹⁰ m²/s, with C_o set at unity) in an infinite domain, (a) shows the general case of how concentration changes as a function of distance in the semi-infinite domain, with time shown as sequential curves. (b) The general case of concentration of a finite nutrient diffusing in space and time, as described by Eq. 21, and as with all 3D graphs herein, depth is shown on the rightward axis and time on the leftward axis. The concentration intensity profile from (b) is shown as a contour plot in distance and time for either (c) D = 10⁻¹⁰ m²/s or (d) D = 10⁻⁹ m²/s, all other parameters being equal. Time is shown until 4000 s in (b–d). (e, f) The condition where, rather than finite diffusant, external C_o remains constant (like ambient oxygen), as described in the semi-infinite domain of Eq. 18. (g, h) Diffusion and concentration profile in a 4-mm-thick finite construct with unlimited supply of oxygen as described in Eq. 22. (e–h) D = 10⁻⁹ m²/s, C_o = 0.22 mM, and time is shown until 10,000 s.

The diffusion equation can also be solved over the semi-infinite distance for diffusant substances that are not consumed by the cells and only have a finite supply in the media (like certain signaling molecules or growth factors). At t = 0, the concentration of diffusant from media into hydrogel appears as a step function at x = 0, where concentration in the media is C_o and concentration in the hydrogel is zero. The solution can be found by treating the finite amount of diffusant at each position as a source of diffusion to adjacent positions, and letting N represent the total number of diffusant particles in the source solution, which in the 1D case is ⁴² \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}N \,=\, \int {C_o} \ d \xi \tag{19}\end{align*} \end{document}

In this system, the variable x is treated as the position where concentration is measured and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\xi$$ \end{document} is the position acting as source of diffusant. The solution over the semi-infinite domain is then the integration over all sources along \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\xi$$ \end{document} , \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}C ( { x , t } ) = \frac { N } { \sqrt { \pi Dt } } \int \nolimits_0^ { \infty } e^ { \left[ \frac { - ( { \xi - x } ) ^2 } { { 4Dt } } \right] } \tag { 20 } \end{align*} \end{document}

producing a solution for C in terms of x and t for the finite nutrient over semi-infinite distance ⁴² (results shown in Fig. 3a–d): \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}C ( { x , t } )\, =\, \frac { N } { \sqrt { \pi Dt } } e^ { - \frac { x^2 } { { 4Dt } } } \tag { { 21 } } \end{align*} \end{document}

Diffusion models of finite slab constructs

Interestingly, similar transformations and solutions using the error function can be used to represent diffusion into a tissue construct of any thickness in the finite domain. Solutions for general diffusion along the x-axis, into a slab, have been derived for various applications,^41,45 which may be modified, in this study, to apply to the particular case of tissue constructs of thickness T: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \begin{split} C ( {x , t} ) &= C_{o} \left[{ \rm erfc} \left(\frac{x} {\sqrt{4Dt}} \right)\right.\\ &\quad - \mathop \sum \limits_{n = 1}^{\infty} ( { - 1} )^{n + 1}{ \rm{erfc}}\left(\frac{({2n}) T - x}{ \sqrt {4Dt}}\right) \\ &\quad \left. + \sum\limits_{n = 1}^{ \infty} ({ - 1})^{n + 1}{\rm {erfc}} \left( \frac{( {2n} ) T + x}{\sqrt {4Dt}} \right) \right] \end{split}\tag{22} \end{align*} \end{document}

In this model, concentration continues to rise even at the deepest point within the construct as long as the nutrient supply remains higher at the source point, as can be seen in the case of oxygen diffusion into a 4-mm-thick slab as graphed in Figure 3g and h. However, this solution still does not completely describe diffusion limitations or consumption of nutrients by cells within the tissue construct.

Steady-state diffusion profiles in metabolically maximized tissue constructs

By substituting the formulas for maximal diffusion depth (Eqs. 4, 10, or 13) into the metabolic concentration profiles for each form of construct (Eqs. 3, 9, or 12), general steady-state solutions for metabolically maximized constructs can be obtained and written in a simplified form: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} { \rm Slab } : C ( x ) = { \frac { \varphi x^2 } { 2D } } - \sqrt { \frac { 2 { C_o } \varphi } { D } } x + { C_o } \tag { { 23\rm a } } \end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} { \rm Cylinder } : C ( r ) = { \frac {\varphi r^2 } { 4D } } \tag { { 23\rm b } } \end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} { \rm Sphere } : C ( r ) = { \frac {\varphi r^2 } { 6D } } \tag { { 23\rm c } } \end{align*} \end{document}

These formulas describe the steady-state metabolic profile through the construct, and it may be verified that these solutions meet the assigned conditions of C = C_o at x = 0 and C = 0 at x = T_max for the slab and C = 0 at r = 0 and C = C_o at r = R_max for the two radial cases. These solutions, however, do not describe the transient phase of unsteady state where concentration profile through the construct progresses toward the steady state from some initial value. The initial concentration of oxygen and nutrients in a newly formed tissue construct is typically zero (C_i = 0), such as when a hydrogel is mixed and cellularized outside of media or under anoxic conditions before cellularization, as is often done with many types of polymer hydrogels, when scaffolds are seeded and transferred to new media, or when a new nutrient is introduced during stages of induction and differentiation. It has been noted that a large amount of cell death occurs in the early initial period immediately after the seeding of a new 3D tissue construct, ⁷⁵ which is likely due to these nutrient deficiencies, and suggests that an understanding of how concentrations change within a tissue construct during transient perturbations is especially important in tissue engineering.

Steady- and nonsteady-state model of slab constructs constrained to maximization with unlimited nutrient

To derive nonsteady-state solutions for constructs with diffusion and metabolism, it is instructive to first derive solutions for the simplest condition for a theoretical maximized construct, where the surface or starting point concentration is held constant (e.g., C = C_o) and the maximal depth concentration is held constant (e.g., C = 0), as would be expected in a maximized construct where all nutrient was completely consumed by the deepest point at steady state. Boundary conditions with constant values are generally known as Dirichlet boundary conditions, and this form of solution is convenient because it provides a foundation for how more complicated solutions will be derived and also because there are many scenarios where the value of each parameter in a construct is not known, meaning that maximal theoretical depth is not known, but because a maximal depth for a viable tissue construct can be observed empirically, this value can simply be inserted in lieu of the individual parameters that affect such limits. Thus, this unique approach can still allow an estimation of diffusion dynamics within the tissue construct even when parameters of diffusivity, initial concentrations of diffusant, cellular density, or metabolic rates of cells are not known. Maximized tissue constructs are also generally of most interest in research applications since they represent the limits of tissue size, metabolic activity, and cell viability, and this method can also allow individual parameters to be interpolated from experimental data or allow modeling of maximal diffusion depth that changes over time or areas of nutrient surplus and deficiency. Such forms of solutions are also useful where consumption of nutrient is maintained by external factors other than cell consumption through the construct, such as when a nutrient source is supplied at one side of a tissue construct and washed out at an opposing side. Such a scenario may occur, for example, when cells are fed with one nutrient media over the tissue surface, but then perfused with a separate media lacking the nutrient underneath the cells to create a gradient through the tissue, or when the nutrient undergoes degradation, sequestration, chelation, absorption, or reaction when exposed at one side of the tissue construct. Solutions where the nutrient molecule is consumed by cells throughout the construct will then be presented in more detail in subsequent sections.

Complete solutions over a finite domain may be obtained by means of a method of separation of variables with the application of specific boundary conditions, a method that was developed for this form of equation by Joseph Fourier in the early 19th century. ⁷⁶ This method assumes that the solution for C(x,t) equals the product of a function of x and a function of t, defined as F(x) and G(t), respectively. \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}C ( {x , t} ) = F ( x ) \;G ( t ) \tag{{24}}\end{align*} \end{document}

The two functions F(x) and G(t) can be plugged into the diffusion Eq. (1) to give \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \frac { 1 } { { DG } } { \frac { dG } { dt } } = \frac { 1 } { F } \frac { d^2F } { dx^2 } \tag { { 25 } } \end{align*} \end{document}

Because the two different functions on the left and right sides of the equation are equal, a solution can be found by setting each side of the equation equal to an arbitrary constant γ and creating a set of ODEs that can be solved with the appropriate boundary conditions over the finite domain: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}\begin{cases} { \frac { dG } { dt } } = \gamma DG \\ \frac { { d^2 } F } { d { x^2 } } = \gamma F \end{cases} \tag { { 26 } } \end{align*} \end{document}

The critical condition of cells consuming exactly enough nutrient that concentration approaches zero at the deepest point in the construct is of particular interest since it describes a maximal tissue size and allows Dirichlet boundary conditions to conveniently be set at both the surface and depth of the tissue construct, as in C(T, t) = 0. The case of a constant supply of oxygen or nutrient at the construct interface is first considered in this study, meaning C(0, t) = C_o. As mentioned previously, the concentration of oxygen and nutrients in the construct at initiation is often zero due to the fact that many tissue constructs are created by seeding cells into hydrogels or scaffolds, meaning C(x, 0) = 0 for the range 0 ≤ x ≤ T at t = 0. If γ > 0 (which can be defined by γ = ɛ ² , where ɛ is the range of eigenvalues \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$ \frac { n \pi } { T } $$ \end{document} , sometimes called Helmholtz eigenvalues after the physicist and physician Hermann von Helmholtz who investigated the form of equation for F above), then \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$F ( x ) = ae^{ \varepsilon x} \,+ \, be^{ - \varepsilon x}$$ \end{document} , which under the boundary conditions above simply results in the uninteresting case of a = b = 0. If γ = 0, then F(x) = ax + b, which under the boundary conditions results in F(0) = C_o = b and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$a \, = \, - \frac { C_o } { T } $$ \end{document} , meaning \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$F ( x ) = C_o ( 1 \, - \, \frac { x } { T } )$$ \end{document} , and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$ \frac { dG } { dt } \,= \, 0$$ \end{document} . This provides a steady-state solution, but because G(t) is a constant in this case, this solution alone does not adequately describe the change in concentration over time. The nonsteady-state solution is found when γ < 0, and the general form of the solutions become \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$F ( x ) = a \cos \left( \frac { n \pi } { T } x \right) + b \sin ( \frac { n \pi } { T } x )$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$G ( t ) = ge^ { - \left( \frac { n \pi } { T } \right) ^2 Dt } $$ \end{document} where a, b, and g are found from the conditions of the system: F(0) = 0 = a, F(T) = b sin (ɛT), meaning that \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$F ( x ) \,= \, \sin ( \frac { n \pi } { T } x )$$ \end{document} . By the principle of superposition then, the linear combination of solutions of F(x) and G(t) can then be given by Fourier series representation of the eigenfunctions: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}C ( x , t ) = C_o \left( 1 - \frac { x } { T } \right) + \mathop \sum \limits_ { n = 1 } ^ { \infty } ge^ { - \left( \frac { n \pi } { T } \right) ^2 Dt } \sin \left( \frac { n \pi } { T } x \right) \tag { 27 } \end{align*} \end{document}

By applying the boundary condition \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$C ( x , 0 ) = 0$$ \end{document} , again meaning that at t = 0 no diffusion has yet occurred, the value of the odd function g can be found by using the orthonormality property of eigenfunctions and multiplying C(x) by F(x) and integrating over the range of x, or by using the related solution given in standard Fourier series tables ⁷⁷ : \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}g = \frac { 2 } { T } \int_0^T - { C_o } \left( { 1 - \frac { x } { T } } \right) \sin \left( \frac { { n \pi } } { T } x \right) dx = { \frac {- 2 C_o } { n \pi } } \tag { { 28 } }\end{align*} \end{document}

It may be noted that the Biot coefficient (B) can be described by the formula \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$B = \frac { hT } { D } $$ \end{document} , where h and D are given in the condition \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$- D \frac { dC } { dx } \,= \, hC$$ \end{document} at the construct surface x = T, and a high Biot coefficient describes a system where the mass transfer coefficient (h) at the surface of the construct is much larger than the diffusivity (D), meaning that a molecule is inhibited at a phasic interface only by diffusivity of the medium and there is not a prescribed time lag for crossing at the interface, and also indicating that the system is defined by Dirichlet-type boundary conditions rather than the Neumann type. Applying the boundary conditions of (1) C(x = 0, t) = C_o, (2) C(x = T, t) = 0, and (3) C(x, t = 0) = 0, a final complete solution for C(x, t) in a maximized finite construct of thickness T thus becomes \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}C ( { x , t } ) = { C_o } \left( { 1 - \frac { x } { T } } \right) - \frac { 2 { C_o } } { \pi } \mathop \sum \limits_ { n = 1 } ^ { \infty } \frac { 1 } { n } e^ { - \left( \frac { n \pi } { T } \right) ^2 Dt } \sin \left( \frac { n \pi } { T } x \right) \tag { { 29 } } \end{align*} \end{document}

The maximal tissue construct thickness obtained from Eq. 4 can be inserted for T, which produces a complete closed-form solution for concentration within a construct in space and time. This solution was graphed in MATLAB as shown in Figure 4a and b, and comparison to numerical PDE solutions in MATLAB demonstrated exact matches for all equivalent conditions (see Appendix for PDE solver methods). The series can be approximated well for early times with n values of about 500 or less, and later times can be adequately represented with an n of only 10 or less. For conditions where the initial concentration in the construct (C_i) is a value other than zero, (C_o − C_i) may be substituted for C_o values and (C − C_i) substituted for C as appropriate. For a model where C_i = 0 and maximal thickness is calculated to be 4 mm, the concentration of oxygen at the halfway point through the construct reaches 50% of its steady-state value in ∼25 min and 90% of its steady-state value in ∼1.2 h.

OPEN IN VIEWER

FIG. 4.

(a, b) Concentration profiles for oxygen in culture media using the constrained solution for unlimited nutrient (Eq. 29) over a construct's spatial domain: C_o = 0.22 mM, D = 10⁻⁹ m²/s, T_max = 4 mm, and time shown from 0 to 4000 s. The view of all axes is shown in (a) and the concentration intensity profile in time and distance (horizontal and vertical axes, respectively) is shown in (b). In (c, d) the concentration profile of glucose from the solution in Eq. 36 is shown using C_o = 10 mM, D = 10⁻¹⁰ m²/s, and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\varphi = 1.25 \times 10^ { - 7 } \frac { mol } { L \cdot s } $$ \end{document} for T_max = 4 mm, until total consumption at t = 3.33 × 10⁵ s at x = 0 with a media volume five times that of the construct volume. Under these conditions, the ramp up in concentration due to diffusion occurs much more quickly than the consumption of nutrient. In (e), the transient ramp up in oxygen concentration within the construct is shown for the first 4000 s as it approaches parabolic quasi-steady state, as described by the solution in Eq. 37 using C_o = 0.22 mM, D = 10⁻⁹ m²/s, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\varphi_ { O_2 } = 2.75 \times 10^ { - 8 } \frac { mol } { L \cdot s } $$ \end{document} , and T_max = 4 mm. The concentration profile of oxygen in the construct over time is shown in (f). The numerical solution with the same conditions as (e, f) is shown in (g, h), as explained in the text. The concentration profile of glucose in the construct from the metabolic solution in Eq. 38 is shown in (i, j) over the valid domain of x = 0 to T_max and until total consumption at 3.3 × 10⁵ s using C_o = 10 mM, D = 10⁻¹⁰ m²/s, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\varphi = 1.25 \times 10^ { - 7 } \frac { mol } { L \cdot s } $$ \end{document} , and T_max = 4 mm, with a media volume five times that of the construct volume.

OPEN IN VIEWER

Steady- and nonsteady-state model of slab constructs constrained to maximization with limited nutrient

For limited nutrient that is consumed by cells from the media (e.g., glucose), it is also important to account for the drop in concentration that occurs in the surrounding media as nutrient is consumed, which in turn feeds back onto the driving source of diffusion. Accurately modeling this behavior with closed-form solutions can be difficult, but the drop in concentration from the media is simply a function of two components: (1) the amount of nutrient consumed by the construct and (2) the amount of nutrient diffused into the construct. By using the total consumption from the tissue construct and adjusting for the difference in volume between the construct and the media volume, C_o can be replaced with a function that describes the rate of metabolic consumption from the media by the construct, where V_c represents the volume of the tissue construct and V_m is the media volume around the tissue construct: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}C_ { media } = { C_o } - \varphi t \frac { V_c } { V_m } \tag { { 30 } } \end{align*} \end{document}

Second, the drop in concentration from the media due to diffusion alone into the construct at any point in time should also be accounted for, represented by \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${ \bar C}$$ \end{document} , with formulas for different construct types given in Eq. 49a–c (although, if desired, for practical simplification this term may also be neglected). Neglecting metabolism for a moment, the molar amount of nutrient in a closed diffusion system remains constant, meaning that \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}C_o V_m = C_{media} {V_m} + { \bar C} {V_c} \tag{{31}}\end{align*} \end{document}

Solving for C_media, the drop in concentration in media due to diffusion alone is described by \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}C_ { media } = { C_o } - { \bar C } \frac { V_c } { V_m } \tag { { 32 } } \end{align*} \end{document}

The formulas for \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${ \bar C}$$ \end{document} in Eq. 49a–c, however, are for a constant unlimited nutrient, and when a limited nutrient diffuses from media into a construct, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${ \bar C}$$ \end{document} cannot rise all the way to C_o since the volume V_m + V_c is greater than V_m alone. Rather, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${ \bar C}$$ \end{document} reaches a new equilibrium in both media and construct, which can be expressed by replacing C_o in Eq. 49a–c with \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$ \frac { C_o V_m } { V_m + V_c } $$ \end{document} , and insertion of this expression along with Eq. 49a–c for \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${ \bar C}$$ \end{document} in Eq. 32 produces a formula for the drop in concentration from media due to diffusion that accounts for the total volume of V_m + V_c: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}C_ { media } = { C_o } - { \bar C } \frac { V_c } { V_m + { V_c } } \tag { { 33 } } \end{align*} \end{document}

Finally, a complete solution for decline of limited nutrient from media due to diffusion and metabolic consumption may then be formed by combining the metabolic consumption and the diffusional compartment shift decrease, replacing C_o in Eq. 29 with \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}C_ { media } = \left( { C_o } - \varphi t \frac { V_c } { V_m } - { \bar C } \frac { V_c } { ( { V_m } + { V_c } ) } \right) \tag { { 34 } } \end{align*} \end{document}

It is important to note that in the case of limited nutrient where C_o is consumed in time, the maximized value of T given by Eq. 4 also becomes a function of time, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}T_ { max_t } = \sqrt { \frac { 2sD \left( C_o - \varphi t \frac { V_c } { V_m } - { \bar C } \frac { V_c } { ( { V_m } + { V_c } ) } \right) } { \varphi } } \tag { { 35 } } \end{align*} \end{document}

which altogether produces the following model for concentration in maximized constructs with limited nutrient:

OPEN IN VIEWER

This solution is valid within the domain of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$T_{{max}_t}$$ \end{document} and until C_o reaches zero in media (given by Eq. 51). The solution shows that the concentration will have a maximal peak at a slightly different time at each point in the construct. Using C_o = 10 mM, D = 10⁻¹⁰ m²/s, and φ for glucose consumption of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$1.25 \times 10^ { - 7 } \; \frac { mol } { L \cdot s } $$ \end{document} , the local maximum at 50% distance into the 4-mm maximal thickness of the construct occurs at 11.6 h with a peak concentration 35% of initial C_o, while at 90% depth into the construct, the local maximum occurs at 8.5 h with a peak concentration 4% of initial C_o. Assuming media volume is only five times that of the construct volume, glucose in the media is completely consumed by 93 h (Fig. 4c, d). At early times, the zero-order consumption plus diffusional equilibration cause an exponential-like decay in nutrient concentration at the construct surface, while at later times once the majority of diffusion has occurred, the decline becomes linear as would be expected with zero-order consumption. The maximal viable thickness of the construct changes as time progresses since the outer concentration of nutrient is changing with time, meaning that nutrient deficiency will occur both in transient states and as nutrient is consumed. The T_max can also be modified to account for changes in φ, such as when cells proliferate and change cell density in the tissue construct.

Steady- and nonsteady-state model of metabolically maximized slab constructs with unlimited nutrient

The prior solutions, however, do not account for the effects of metabolism by cells homogenously distributed through a construct. For any case where metabolism of nutrient by cells is greater than zero (φ > 0), the concentration profile for the nutrient through the construct will be parabolic, as described previously. A complete steady- and nonsteady-state model for diffusion and metabolism with unlimited nutrient may be derived as before by again implementing superposition of particular solutions that satisfy the boundary conditions of (1) C(x = 0, t) = C_o, (2) C(x = T_max, t) = 0, and the initial state of (3) C(x, t = 0) = 0, and satisfy the known steady-state solution (Eq. 3): \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \begin{split}C ( x , t ) = &C_o + { \varphi x^2 \over 2D} - { \varphi Tx \over D} - {2C_{o} \over \pi} \mathop \sum \limits_{n = 1}^{ \infty} {1 \over n} e^{ - \left( n \pi \over T \right) ^2 Dt} \\&\sin \left( n \pi {2 \varphi Tx - \varphi x^2 \over 2C_oD} \right)\end{split} \tag{37}\end{align*} \end{document}

The T_max value, either calculated from parameters or observed empirically, can then be inserted for T. A graph of Eq. 37 with typical oxygen parameters is shown in Figure 4e and f. Comparison of Eq. 37 to numerical solutions under the same conditions shows that the analytic and numerical solutions are identical at steady state, but there is a small variance at early transient times in the deepest point in the construct, where concentration briefly becomes negative in numerical models due to the construct having a constant φ, but nutrient not having yet reached the full depth of the construct, as demonstrated in Figure 4g and h. This transient negative deficiency is only present in the numerical models due to the assumed condition that consumption of nutrient continues even in the absence of nutrient, thereby building up a deficiency until nutrient arrives, as opposed to the analytic models which keep nutrient concentration at the initial zero value until nutrient has diffused to the space, as would be measured in a culture system. Yet the analytic models presented herein still closely replicate the numerical models with all appropriate boundary conditions and may in fact represent more realistic conditions where cells cannot metabolize nutrient until the nutrient arrives, where nutrient concentration cannot be measured as a negative value in culture, and where cells do not simply restore nutrient deprivation as a linear summation over time. Mathematical code for both the numerical and analytical approaches is provided for convenience in the Appendix.

Steady- and nonsteady-state model of metabolically maximized slab constructs with limited nutrient

As before, for limited nutrient that is consumed by cells from the media, C_o can be replaced with Eq. 34, where \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${ \bar C}$$ \end{document} is given by Eq. 49a and the maximized value of T becomes a function of time (Eq. 35), forming a complete model for concentration within the construct of limited nutrient that is consumed from the media: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \begin{split}C ( x , t ) =& \left( C_o - \varphi t {V_c \over V_m} - { \bar C} {V_c \over \left( V_m + V_c \right) } \right) + { \varphi x^2 \over 2D} \\&- x \sqrt {2 \varphi \left( C_o - \varphi t {V_c \over V_m} - { \bar C} {V_c \over \left( V_m + V_c \right) } \right) \over D} \\&- {2 \left( C_o - \varphi t {V_c \over V_m} - { \bar C} {V_c \over \left( V_m + V_c \right) } \right) \over \pi} \mathop \sum \limits_{n = 1}^{ \infty} {1 \over n} e^{ - \left( {n \pi \over T} \right) ^2 Dt} \\&\sin \left( n \pi x { \sqrt{8D \varphi \left( C_o - \varphi t {V_c \over V_m} - { \bar C} {V_c \over \left( V_m + V_c \right) } \right) } - \varphi x \over 2D \left( C_o - \varphi t {V_c \over V_m} - { \bar C} {V_c \over \left( V_m + V_c \right) } \right) } \right)\end{split} \tag{38}\end{align*} \end{document}

The result of Eq. 38 with typical glucose parameters is shown in Figure 4i and j.

Diffusion models of finite cylindrical constructs

By the same method of separation of variables used to derive the 1D solution, the general equation can be obtained for the cylindrical case, where again F is a function of radial distance and G is a function of time: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \frac { 1 } { DG } \frac { dG } { dt } = \frac { 1 } { rF } \frac { d } { dr } \left( r \frac { dF } { dr } \right) \tag { 39 } \end{align*} \end{document}

The general solutions to the form of F(r) in this case involve Bessel functions, and with the boundary conditions of constant concentration of nutrient in the media, C(R, t) = C_o, and initial concentration of zero within the cylinder, C(r, 0) = 0 for the range 0 ≤ x ≤ R, as well as the assumption that the length of the cylinder is significantly greater than the radius of the cylinder, the solution is found by previously described methods^39,41,45 \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}C ( r , t ) = C_o \left[ 1 - \frac { 2 } { R } \mathop \sum \limits_ { n = 1 } ^ { \infty } e^ { - ( \lambda_n ) ^2 Dt } \frac { J_o ( r \lambda_n ) } { \lambda_n J_1 ( R \lambda_n ) } \right] \tag { 40 } \end{align*} \end{document}

where J₀ and J₁ are Bessel functions of the first kind (orders zero and one, respectively), and the eigenvalues λ_n are the roots of J₀(Rλ_n) = 0. By using the values in Table 4, where the first 10 roots are shown, λ_n can be solved and inserted in the summation series. Bessel functions can be approximated with a simplified formula ⁷⁸ : \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \begin{split}J_0 ( x ) \approx \sqrt {{2 \over \pi x}} \left( 1 - {1 \over 16{ \rm x}^2} + {53 \over 512{ \rm x}^4} \right) \cos \left( x - { \pi \over 4} - {1 \over 8x} + {25 \over 384x^3} \right)\end{split} \tag{41\rm a}\end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \begin{split}J_1 ( x ) \approx \sqrt {{2 \over \pi x}} \left( 1 + {3 \over 16{ \rm x}^2} - {99 \over 512{ \rm x}^4} \right) \cos \left( x - {3 \pi \over 4} - {3 \over 8x} + {21 \over 128x^3} \right)\end{split} \tag{41\rm b}\end{align*} \end{document}

Therefore, although quite unwieldy, an approximate closed-form solution for cylindrical diffusion can be written as follows, where values for λ_n are given in Table 4:

OPEN IN VIEWER

The primary weakness in this approach, however, is that the approximations of Bessel functions become inaccurate at small values of r, and the sum function is inaccurate for small values of t unless carried out to large n values, meaning that at early times in the center of the sphere, this simplified closed-form solution falls apart. However, by using actual Bessel function calculations in MATLAB, more accurate models for general cylindrical diffusion can be achieved (Fig. 5a–d).

OPEN IN VIEWER

FIG. 5.

(a–d) Cylindrical diffusion profiles from Eq. 40 along the radial axis through time assuming no metabolism. The Bessel function solutions are shown in a cylindrical construct of 1 mm radius for the case of unlimited oxygen diffusion (a, b) and glucose diffusion (c, d); parameters for oxygen were C_o = 0.22 mM and D = 10⁻⁹ m²/s with time shown until 1000 s and for glucose were C_o = 10 mM and D = 10⁻¹⁰ m²/s. Cylindrical diffusion profiles of oxygen diffusion in a construct with a maximized diameter of 2 mm from Eq. 43 are shown in (e, f); parameters were C_o = 0.22 mM, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\varphi_ { O_2 } = 8.8 \times 10^ { - 7 } \frac { mol } { L \cdot s } $$ \end{document} , and D = 10⁻⁹ m²/s. The numerical solution with the same conditions as (e, f) is shown in (g, h), as explained in the text. Cylindrical diffusion profiles of limited nutrient (glucose) consumed in a construct with a maximized diameter of 2 mm described by Eq. 44 over the valid domain of r = 0 to R_max and until C_o reaches zero in media are shown in (i, j) assuming the media volume is only five times more than the cylindrical construct volume, with standard glucose parameters of C_o = 10 mM, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\varphi_ { Gluc } = 4.0 \times 10^ { - 6 } \frac { mol } { L \cdot s } $$ \end{document} , and D = 10⁻¹⁰ m²/s.

OPEN IN VIEWER

Steady- and nonsteady-state model of metabolically maximized cylindrical constructs with unlimited nutrient

As with the 1D case, the 2D cylindrical case also uses the steady state defined by metabolism in the construct (Eq. 12). The cylindrical equation for a maximized value of R that results in zero concentration at the center of the cylinder is implemented with boundary conditions of (1) C(r = R_max, t) = C_o, (2) C(r = 0, t) = 0, and the initial state of (3) C(r, t = 0) = 0, and a complete steady- and nonsteady-state model becomes \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}C ( r , t ) = \frac { C_o r^2 } { R^2 } + \frac { 2C_o } { \pi } \left[ \mathop \sum \limits_ { n = 1 } ^ { \infty } \frac { - 1^n } { n } e^ { - \left( \frac { n \pi } { R } \right) ^2 Dt } \sin \left( n \pi \frac { r^2 } { R^2 } \right) \right] \tag { 43 } \end{align*} \end{document}

Graphical representation of this equation is shown in Figure 5e and f. It is useful to note that in the radial cases, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$ \frac { C_o r^2 } { R_ { max } ^2 } $$ \end{document} equals \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$ \frac { \varphi r^2 } { 4D } $$ \end{document} for the cylinder (Eq. 23b) or \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$ \frac { \varphi r^2 } { 6D } $$ \end{document} for the sphere (Eq. 23c). If desired, the concentration profile through the construct may be compared to the 1D slab case by reversing the spatial axis by substitution of (R_max − r) in place of r. The numerical solution under the same conditions of Eq. 47 and with the same parameters used in Figure 5e and f is shown in Figure 5g and h, demonstrating a transient negative nutrient concentration at an early phase where cells are presumed to continue nutrient consumption even as nutrient diffuses toward the inner region of the construct, as described in Steady- and Nonsteady-State Model of Metabolically Maximized Slab Constructs with Unlimited Nutrient section.

Steady- and nonsteady-state model of metabolically-maximized cylindrical constructs with limited nutrient

A model for limited nutrient consumed in the maximized construct can then be created as before by substituting \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$C_o = \left( C_o - \varphi t \frac { V_c } { V_m } - { \bar C } \frac { V_c } { \left( V_m + V_c \right) } \right)$$ \end{document} from Eq. 34 and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$R_ { max_t } = \sqrt \frac { 4D \left( C_o - \varphi t \frac { V_c } { V_m } - { \bar C } \frac { V_c } { \left( V_m + V_c \right) } \right) } { \varphi } $$ \end{document} from Eq. 35 into Eq. 43 (where \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${ \bar C}$$ \end{document} is given by Eq. 49b and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$R_{max_t}$$ \end{document} is substituted in place of R), which simplifies substantially after canceling terms. Because these models describe tissue constructs at a constant thickness or constant radius, and because the maximal depth that changes over time should be measured from the surface of the constructs (as was shown in the 1D case of a slab with limited nutrient), the variable r in radial cases of limited nutrient should also be replaced with \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\left( r - R_{max} + R_{max_t} \right)$$ \end{document} as shown in Eq. 44.

A graph of Eq. 44 over the valid domain of r = 0 to \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$R_{max_t}$$ \end{document} and until C_o reaches zero in media is shown in Figure 5i and j, where the V_m was again assumed to be five times that of V_c. Such a model is relevant not only for tissue constructs but also for bioreactor systems where cells are held in cylindrical culture constructs and perfused with a concentration gradient of nutrient.

Diffusion models of finite spherical constructs

Spherical diffusion solutions can be formed from the 3D diffusion equation by using the separation of variables method as shown in the linear and cylindrical cases. When independent functions of time and distance are substituted into the spherical diffusion equation, the following is obtained: OPEN IN VIEWER \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \frac { 1 } { DG } \frac { dG } { dt } = \frac { 1 } { r^2 F } \frac { d } { dr } \left( r^2 \frac { dF } { dr } \right) \tag { 45 } \end{align*} \end{document}

The general solutions produced are then \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$F ( x ) = \frac { a } { r } \cos ( \varepsilon r ) + \frac { b } { r } \sin ( \varepsilon r )$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$G ( t ) = ge^{ - \varepsilon^2 Dt}$$ \end{document} . The condition of C(r = 0, t) = 0 is satisfied by letting a = 0, meaning \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$F ( x ) = b\, \sin ( \varepsilon r )$$ \end{document} . Using the boundary conditions of constant nutrient at the surface of the sphere, finite nutrient levels at the center of the sphere, and no initial nutrient within the sphere as described previously, C(R, t) = C_o and C(r, 0) = 0 for the range 0 ≤ x ≤ R, a complete diffusion solution for a general spherical construct without metabolic consumption can be written as ⁴¹ \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}C ( r , t ) = C_o + \frac { 2RC_o } { \pi r } \left[ \mathop \sum \limits_ { n = 1 } ^ { \infty } \frac { - 1^n } { n } e^ { - \left( \frac { n \pi } { R } \right) ^2 Dt } \sin \left( \frac { n \pi r } { R } \right) \right] \tag { 46 } \end{align*} \end{document}

This formula shows that in a nonmaximized 2-mm-diameter sphere with no consumption of oxygen and D = 10⁻⁹ m²/s, it takes just less than 10 min for diffusion of oxygen at the center of the sphere to reach 99% of surrounding oxygen levels (Fig. 6a–d). Likewise, for a constant ambient glucose with negligible consumption of nutrient from media (e.g., large V_m:V_c ratio) and D = 10⁻¹⁰ m²/s, glucose reaches 99% of surface levels at the center of a maximized sphere in just under 1.5 h.

OPEN IN VIEWER

FIG. 6.

 (a–d) Concentration profile for spherical diffusion from Eq. 46 of oxygen (a, b) using R = 1 mm, C_o = 0.22 mM, D = 10⁻⁹ m²/s, and assuming no metabolism to 1000 s and glucose (c, d) using R = 1 mm, C_o = 10 mM, D = 10⁻¹⁰ m²/s to 4000 s. The left axis represents time and the right axis represents the radial distance into the sphere, where zero is the center of the sphere. Thus, if the right axis alone was rotated about its origin in all spatial directions, it would represent the dimensions of the entire sphere. In (e, f) concentration profiles are shown for unlimited nutrient (oxygen) from Eq. 47 for a 2 mm diameter metabolically maximized sphere ( \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\varphi_ { O_2 } \, = \, 1.32 \times 10^ { - 6 } \frac { mol } { L \cdot s } $$ \end{document} ) where nutrient is fully consumed by the center of the construct, while remaining constant at the surface of the sphere. The numerical solution with the same conditions as (e, f) is shown in (g, h), as explained in the text. In (i, j) concentration profiles of a 2 mm maximal diameter sphere with limited nutrient (glucose) and zero-order metabolism described by Eq. 48 over the valid domain of r = 0 to R_max and until C_o reaches zero in media using C_o = 10 mM, D = 10⁻¹⁰ m²/s, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\varphi_ { Gluc } \, = \, 6.0 \times 10^ { - 6 } \frac { mol } { L \cdot s } $$ \end{document} , and assuming the media volume is only five times greater than the spherical construct volume.

OPEN IN VIEWER

Steady- and nonsteady-state model of metabolically maximized spherical constructs with unlimited nutrient

As described previously, the steady-state concentration is determined by metabolism and diffusivity in the construct. The spherical equation for a maximized value of R that results in zero concentration at the center of the cylinder is again implemented with boundary conditions of (1) C(r = R_max, t) = C_o, (2) C(r = 0, t) = 0, and (3) C(r, t = 0) = 0, and a complete steady- and nonsteady-state model then becomes \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}C ( r , t ) = \frac { C_o r^2 } { R^2 } + \frac { 2C_o } { \pi } \left[ \mathop \sum \limits_ { n = 1 } ^ { \infty } \frac { - 1^n } { n } e^ { - \left( \frac { n \pi } { R } \right) ^2 Dt } \sin \left( n \pi \frac { r^2 } { R^2 } \right) \right] \tag { 47 } \end{align*} \end{document}

As before, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\frac{C_o r^2}{R_{max}{^2}}$$ \end{document} simplifies to \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$ \frac { \varphi r^2 } { 6D } $$ \end{document} , as given in Eq. 23c, and, if desired, the concentration profile may be compared to the 1D case by substituting (R_max − r) for r. The results of Eq. 47 are shown in Figure 6e and f for a 2 mm maximal diameter sphere where drops in nutrient at the surface of the sphere are negligible, but nutrient is fully consumed by the center of the construct. The results for this spatial concentration profile are also supported by experimental evidence for diffusion in spherical constructs.^72,79–81 The numerical solution under the same conditions of Eq. 47 and with the same parameters used in Figure 6e and f is shown in Figure 6g and h, demonstrating a transient negative nutrient concentration at an early phase where cells are presumed to continue nutrient consumption even as nutrient diffuses toward the inner region of the construct, as described in Steady- and Nonsteady-State Model of Metabolically Maximized Slab Constructs with Unlimited Nutrient section.

Steady- and nonsteady-state model of metabolically maximized spherical constructs with limited nutrient

A model for limited nutrient consumed in the maximized construct can then be created as before by substituting \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$C_o = \left( C_o - \varphi t \frac { V_c } { V_m } - { \bar C } \frac { V_c } { \left( V_m + V_c \right) } \right)$$ \end{document} from Eq. 34 and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$R_ { max_t } = \sqrt \frac { 6D \left( C_o - \varphi t \frac { V_c } { V_m } - { \bar C } \frac { V_c } { \left( V_m + V_c \right) } \right) } { \varphi } $$ \end{document} from Eq. 35 into Eq. 47 (where \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${ \bar C}$$ \end{document} is given by Eq. 49c and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$R_{max_t}$$ \end{document} is substituted in place of R). As before, because these models describe constructs at a constant radius and because the maximal depth that changes over time should be measured from the surface of the constructs, the variable r in radial cases of limited nutrient should also be replaced with \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\left( r - R_{max} + R_{max_t} \right)$$ \end{document} , as shown in Eq. 48.

OPEN IN VIEWER

A graph of Eq. 48 over the valid domain of r = 0 to \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$R_{max_t}$$ \end{document} and until C_o reaches zero in media is shown in Figure 6i and j. A summary of all models for diffusion and metabolism of limited and unlimited nutrients in slabs, cylinders, and spheres is provided in Table 5.

Models of whole-construct nutrient consumption from media supply

A model of the amount of an unlimited nutrient diffused into a construct can be made by integrating the molar transfer across the outer surface area of a construct from time zero until the time of interest,^41,45,46 as represented by the following formulas showing the material balance diffused into the construct at any point in time assuming no metabolism of nutrient: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} { \rm Slab } : { \bar C } = C_o \left[ 1 - \frac { 8 } { \pi^2 } \sum \nolimits_ { n = 1 } ^ { \infty } \frac { 1 } { ( 2n - 1 ) ^2 } e^ { - \left( \frac { ( 2n - 1 ) \pi } { 2T } \right) ^2 Dt } \right] \tag { 49\rm a } \end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} { \rm Cylinder } : { \bar C } = C_o \left[ 1 - \sum \nolimits_ { n = 1 } ^ { \infty } \frac { 4 } { ( \lambda_n R ) ^2 } e^ { - \lambda^2_n Dt } \right] \tag { 49\rm b } \end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} {\rm Sphere } : {\bar{C}} = C_o \left[ 1 - \frac{6}{\pi^{2}} \sum\nolimits_{n = 1 }^{ \infty } \frac{1}{n^2} e^{-( \frac{ n \pi}{ R} )^2 Dt } \right] %\tag { 49\rm c } \end{align*} \end{document}

where λ_n of the cylindrical equation are given by the roots of J_o(Rλ_n) = 0 listed in Table 4. For limited nutrient in media, the C_o will drop as the nutrient diffuses into the tissue construct, with the new equilibrium concentration becoming \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$ \frac { C_o V_m } { \left( V_m + V_c \right) } $$ \end{document} , which can be substituted for C_o above to determine final equilibrium concentration in the system. The condition of insufficient nutrient supply for a whole tissue construct, as measured from the media, can then be described by accounting for the limited nutrient consumed and diffused from ambient media. Under the conditions of constant φ, the amount of consumption by whole construct metabolism up until any time t is φtV_c. A critical threshold below which the function or viability of construct is impaired or threatened (C_critical) can be assigned to the system, which must typically be found empirically. Assuming quasi-steady-state diffusion is reached before the nutrient is completely metabolized (which is typically the case), the following formula then allows investigation of the time until critical cell dysfunction within the construct: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}C_ { media } = C_o \frac { V_m } { \left( V_m + V_c \right) } - \varphi t \frac { V_c } { V_m } \tag { 50 } \end{align*} \end{document}

where C_media is compared to C_critical. By setting C_media equal to zero, the time to total consumption can also be found: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}t_ { consum } = \frac { C_o V_m^2 } { \varphi ( V_c V_m + V_c^2 ) } \tag { 51 } \end{align*} \end{document}

Also of note, in models of maximized tissue constructs with a critical threshold of nutrient above zero, the maximal construct dimensions may simply be determined for any dimensionality s by \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}R_ { max } = \sqrt \frac { 2sD \left( C_o - C_ { critical } \right) } { \varphi } \tag { 52 } \end{align*} \end{document}

Timing and effects of nutrient consumption

Although the diffusivity of oxygen can be ∼10 times higher compared with glucose, the available concentration of oxygen is typically one to two orders of magnitude less than glucose (glucose is typically 10–20 mM in organoid media) and the molar consumption rate of oxygen is generally higher. Thus, although glucose may take longer to reach steady state in a maximized tissue construct due to a lower diffusion coefficient, the higher concentration can still enable glucose to quickly reach adequate levels to support cell viability. Further complexity is introduced by the fact that the cell viability will depend not only on a critical threshold of nutrient concentration but also on the duration of time that the cell is deprived of the nutrient, by the fact that different cells have different thresholds of viability for different nutrients and one nutrient can also influence availability of other nutrients. For example, when a nutrient like glucose is low, other forms of nutrients can also be used as a compensatory energy supply, or if oxygen levels are low, anaerobic mechanisms can compensate with certain limitations. The Crabtree effect describes an increase in oxygen availability during concomitantly high glucose, or a decrease in oxygen availability during low glucose, and insufficient oxygen can likewise cause higher glucose consumption.^82–84 Therefore, the determination of whether the concentration drops below a critical viability threshold must be made based on the characteristics of the cell types and the conditions of the system.

The determination of which nutrient is the limiting factor in a construct can be made by comparing thresholds and the concentration profiles between nutrients. For example, the values of C(r, t) for each nutrient (either steady state or steady and unsteady states for any point in time) can be subtracted from each other or graphed simultaneously and compared to determined critical thresholds. The effective R_max can be calculated from the nutrient that gives the more limiting value, or, when individual parameters are not known, the observed R_max for each type of nutrient can be substituted into the solutions separately for each nutrient to compare the concentration profiles. The models proposed above show glucose to initially be more of a temporary limiting factor for cell metabolism than oxygen due to its lower diffusivity, but after diffusion has approached steady state, oxygen is the limiting nutrient in all cases under the given range of parameters. However, glucose limitations can also become comparable to oxygen limitations when upper limits for glucose and lower limits for oxygen are used (i.e., oxygen consumption and glucose diffusivity are minimized and glucose consumption and oxygen diffusivity are maximized).

The measured cell density in cerebral organoids was found to be \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\rho = 4.82 ( { \pm 1.41} ) \times 10^{11} \ cells / L$$ \end{document} . With this cell density, at the upper range of oxygen diffusivity and the lower range of neural metabolic oxygen consumption, the maximal diameter of organoids is predicted to be 1.4 mm at best, while typical glucose parameters suggest a diameter of about 4 mm, and even stringent glucose parameters suggest a diameter of at least 1.4 mm, thus demonstrating how oxygen is more of a limiting nutrient than glucose at steady state. Under the models presented herein, with R_max dictated by oxygen and glucose consumption set in the normal range, glucose concentrations are maintained at viable levels at all depths and at all times except in the early transient stage after placing cells in hydrogel—a direct comparison of concentrations using the difference between the complete steady- and unsteady-state spherical solutions for glucose and oxygen (with glucose set at 10 mM and oxygen at 0.22 mM) shows that glucose concentration remains lower than oxygen concentration for only the first 8 min. Even though oxygen may not be as limiting as glucose in the early phase of diffusion into a tissue construct, it is plausible that the combination of low oxygen and nutrients at initiation could produce a compounding insult, threatening viability in early states of tissue formation. ⁷⁵ As described, however, oxygen generally becomes the primary limiting factor at steady states in metabolically active constructs.

It is noteworthy that during periods of high neural activity, rapid spiking can increase oxygen consumption by approximately fivefold over baseline firing rates, ⁵¹ and it appears that transient states of glycolysis occur over oxidative metabolism. ⁸⁵ Astroglia compose about half of the number of cells in the brain,^86,87 and their activity has been estimated to account for about 15% of total brain energy consumption. ⁸⁵ Astrocyte metabolism may also increase during the neural activity due to glutamate recycling and astrocytes may take up a significantly higher percentage of glucose and shuttle it as energy to neurons in the form of lactate.^85,88 These models do not account for metabolic coupling or metabolite production; in particular, lactic acid and carbon dioxide are metabolic products that locally increase acidity in solution, which can in turn affect material properties and cell viability, and when oxygen concentration is low or cannot meet cellular demands, lactic acid production can increase through anaerobic glycolysis. In normal neural tissue, a higher neural activity is also coupled to local increases in blood flow, which serves to restore external concentration of surrounding nutrient to initial levels. Thus, it may be important to account for these phenomena within a tissue construct in certain scenarios, but otherwise it can be assumed that nutrient transport systems are independent.

Application of diffusion models to cerebral organoid spheroids

As described in the previous section, the maximal predicted diameter of cerebral organoids (without central cell death) is 1.4 mm. This predicted diameter is based on the physics of diffusion, the observed cell density, and a lower range of reported metabolic consumption rates for oxygen. A low metabolic activity is plausible since the spontaneous neural activity appears to be infrequent in organoids, with a spike activity occurring only every few minutes on average in mature organoids. ³¹ In addition, pluripotent cells are known to sustain quiescent states and may favor nonoxidative glycolysis over oxidative phosphorylation—it has been shown that preimplantation mouse embryos tend to make ATP primarily by oxidative phosphorylation, while the implanted embryos favor glycolysis, and similarly, naive pluripotent stem cells tend to mix oxidative phosphorylation and glycolysis, while more primed pluripotent stem cells lean toward glycolysis ¹⁶ ; other reports have suggested that glycolysis predominates in metabolic profiles of human pluripotent and hematopoietic stem cells, which may not only be a consequence of stem-like states but may also actively help mediate and maintain such states.^{1,2,8,15,19,89} However, metabolic states during stages of neural differentiation are currently not well understood.

Nevertheless, even at the lowest range of reported neural oxygen consumption rates and with known parameters of cell density, diffusivity, and metabolic rate of consumption, the maximal size of the organoid was still able to exceed the predicted maximal size. Cerebral organoids grown from hiPSCs reached an average diameter of 1.8 mm by day 40 (Figs. 7 and 8). Although many of the cerebral organoids obtained widths up to 4 mm along the longest axis, organoids typically had shorter distances between the deepest points and the nearest surface, which shortens the distance needed to perfuse the construct, and it was this shortest distance to the deepest cells that was measured. Notably, the growth characteristics of organoids closely follow those predicted by numerically solved diffusion-limited spherical growth models,^47–49 as shown in Figure 7, meaning that organoid growth is limited by diffusion or by the metabolites produced by diffusion limitations. This, along with the fact that decreased viability is known to occur at the center of organoids, ³¹ supports the idea that diffusion depths are indeed the limiting factor in construct size.

OPEN IN VIEWER

FIG. 7.

Graph showing the average diameter of cerebral organoids over time (± standard deviation). Organoids were able to exceed the maximal diameter predicted by typical oxygen consumption rates in homogenously dense neural constructs, but a model of dual layer cell densities in a spherical construct demonstrates a way of achieving an enhanced maximal diameter. A model of diffusion-limited growth was also applied and fitted to the data, demonstrating that growth follows expected behavior of a diffusion-limited system.

OPEN IN VIEWER

FIG. 8.

Induced pluripotent stem cells were cultured and passaged on a feeder-independent matrigel surface (a). Organoids demonstrated a distinct early outer layer of neuroepithelium that exhibited various forms of architecture (day 20 shown), including broad neuroepithelium (arrowheads in b) or compact localized formations of expanded neuroepithelia resembling 3D rosettes (arrows in c). In (d–f), stained sections of a 40-day old organoid demonstrated a dense outer region with neural identity [(d) green = TUJ1 β-III-tubulin, (e) blue = Hoechst stain, (f) merged], which according to the presented model may serve to enhance organoid size and diffusion characteristics.

There are several explanations for why a discrepancy in predicted versus experimental maximal diameters arose, including the idea that at least a portion of cells may be in a low metabolic state or may be using nonoxidative energy sources, meaning that \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\varphi_{o_2}$$ \end{document} would be lower than predicted from reported neuronal metabolism values. It is also possible that cells may continue to expand the organoid in the outer regions even while cells at the center suffer ischemia. Also, because organoids are not perfect spheres, not all diffusion must occur through the measured density of metabolically active cells, or in other words, the number of cells found in the organoid overrepresents how many cells would be present in a sphere of the same measured radius, and therefore the R_max may have been underestimated by Eq. 10 due to an overestimated effective cell density. Importantly, it is also recognized that the internal structure of organoids is not homogenous, with dense cortical-like layers in outer regions and less dense progenitor regions internally. It was therefore hypothesized that regional architecture within the construct might also influence diffusion characteristics and account for observed discrepancies, and a simple model was created to investigate this question.

Multicompartment model and basis for neural migration and cortical development in avascular spheres

An important question is whether there is an optimal architecture and regional cell density in the organoid that would maximize its potential size. To answer this question, a model was derived to represent how R_max can change when compartmentalizing cells into an inner and outer shell of two different densities given the same amount of total cells in the construct and same nutrient diffusivity. In this model, R_max is first set for a constant homogenous density of cells as described in Eq. 10, with any values for parameters φ, C_o, and D. The sphere is then divided into two concentric spheres, where ϒ represents the proportional distance (from 0 to 1) where the interface between the two shells is located along the radial axis to R_max, and Ω represents the proportion of cells (from 0 to 1) within the volume of the outer shell between ϒR_max and R_max. The densities of cells in the outer and inner shells (called ρ₂ and ρ₁, respectively) can be expressed in terms of the original density (ρ₀), and because it is assumed that cellular metabolism rates stay the same regardless of where cells are located (although regional variations in metabolism and diffusivity could also be integrated into this model), the variable φ₂ can then be used to represent the higher consumption through the new higher density in the outer shell relative to the original φ₀, and the variable φ₁ can be used to represent the lower metabolic consumption within the inner shell relative to φ₀ (assuming, for simplicity, that the shift of cells into a new density in the outer shell occurs before expansion of the sphere and the inner core can then adjust its radius based on the new level of available oxygen): \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}\varphi_2 = \frac { \Omega } { \left( 1 - \Upsilon { ^3 } \right) } \varphi_0 \tag { 53 } \end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}\varphi_1 = \frac { ( 1 - \Omega ) } { \Upsilon ^3 } \varphi_0 \tag { 54 } \end{align*} \end{document}

The drop in concentration across the outer shell can then be found from the following formula, where C_s represents the concentration at the outer rim of the inner sphere, that is, concentration at the interface between the two shells: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}C_s = \frac { \varphi_2 } { 6D } \left( ( \Upsilon R_ { max } ) ^2 - R_ { max } { } ^2 \right) + C_o \tag { 55 } \end{align*} \end{document}

The new potential maximal radius of the inner shell, designated as R_s, can be described with the newly modified density in the inner shell, φ₁, and the new initial concentration at the outer rim of the inner shell, C_s: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}R_s = \sqrt \frac { 6DC_s } { \varphi_1 } \tag { 56 } \end{align*} \end{document}

Finally, the radial depth of the modified outer and inner shells can be summed to find a new optimized maximal radius, designated as R_opt: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}R_{opt} = R_s + ( 1 - \Upsilon ) R_{max} \tag{57}\end{align*} \end{document}

Substituting Eq. 53 into Eq. 55, Eq. 54 into Eq. 56, Eq. 10 into Eqs. 55 and 57, Eq. 55 into Eq. 56, and Eq. 56 into Eq. 57, the following relationship emerges for finding an optimized maximal radius as a function of regional architecture modifications in the sphere: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \frac { R_ { opt } } { R_ { max } } = ( 1 - \Upsilon ) + \sqrt { \left( \frac { \Upsilon ^3 } { 1 - \Omega } \right) \left( \frac { \Omega ( \Upsilon ^2 + 1 ) } { ( 1 - \Upsilon ^3 ) } + 1 \right) } \tag { 58 } \end{align*} \end{document}

This solution is graphed in Figure 9 and shows that when as many of the cells as possible migrate into as thin of outer shell as possible, the sphere can obtain a maximal radius several times greater than if the same number of cells were homogenously distributed throughout the sphere. As an example, if 90% of cells (Ω = 0.9) in a homogenously cellularized sphere are moved into the outer half-volume of the sphere (defined at \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$r = R / \root3\of2$$ \end{document} or ϒ = 0.79), then the new achievable maximum radius (R_opt) is 1.5 times greater than the original R_max. In these cerebral organoids, the observed maximal diameter of 1.8 mm (which is 1.3 times greater than the 1.4 mm limit suggested by the homogenous cellular density model) would suggest that ∼85% of cells (Ω = 0.85) must be in the outer half-volume of the sphere (ϒ = 0.79), as opposed to the 50% that would be expected in the same volume in a homogenous construct, which matches well with the observation that neural precursor cells migrate outwardly in early stages to form high density cortical-like structures in both the brain and organoids^31,90,91 (also see Fig. 8). Clearly this is only a simplified model that neglects the fact that as regions shift and expand, concentration profiles, shell thicknesses, and cell densities will change in response to each other, and mechanical limitations prohibit the formation of an infinitely thin outer shell of cells, but this model demonstrates the concept that diffusion through a thin but dense layer of larger surface area is more efficient than diffusion through a sparse but thick layer of smaller surface area.

OPEN IN VIEWER

FIG. 9.

A multicompartment spherical model for cerebral organoids is shown (a). Curve A demonstrates the general concentration profile for any metabolized gas or nutrient and the maximal diameter of 1.4 mm (dashed line) for a homogenous distribution of cells, while curve B demonstrates the concentration profile and modified maximal radius when a percentage of cells shifts into a denser outer rim (zone 1), all other parameters remaining constant, thereby enabling an enhanced maximal diameter of 1.8 mm (solid border). The percentage of ambient nutrient concentration within the construct is also marked on the concentration profile curves. The maximal diameter is enhanced when more metabolically active cells localize densely into zone 1. Zone 2 may represent an intermediate region of progenitor cells, zone 3 may represent a region of multipotency preserved by hypoxic conditions, and zone 4 represents a potentially ischemic region. (b) A graph of Eq. 58 demonstrating how the maximal achievable radius of a spherical construct increases when a homogenously cellularized sphere is modified into a two-shell sphere with a higher percentage of the cells (Ω) placed into the outer shell beginning at point ϒ along the radius between 0 and R_max.

It is often suggested that a larger and more convoluted cortex confers the advantage of a larger surface area for cortical networks, but it is not explained why neurons must be on the outer surface in the first place. The above model is an algebraic demonstration showing how the largest potential size of a cellularized sphere is achieved by localizing the majority of metabolically active neurons into a thin and dense outer shell, and it may hold many implications for why inchoate and avascular neural systems develop with neural precursors migrating outwardly. A multicompartment model may serve to illustrate internal architecture and subspecialized functions of cells within cerebral organoids, as shown in Figure 9. At the highest oxygen concentrations, metabolically active cells like neurons can best function in a dense outer zone (zone 1) near a high oxygen supply, and this external localization also best enables growth of the organoid and corresponding expansion of neurons in this region. In the second zone, marginally lower oxygen levels may promote expansion of neuroglial precursor populations, a phenomenon which has been observed in cell cultures with controlled hypoxia levels.^11–13 In the central region of still lower oxygen levels, stem cell pluripotency may be maintained as has also been observed in cell cultures with controlled hypoxia levels,^4–6,14 although other evidence has also suggested that hypoxia can induce differentiation.^9,10,18 Finally, zone 4 represents a potential hypoxic region where cell function and viability are threatened.