A triphasic constrained mixture model of engineered tissue formation under in vitro dynamic mechanical conditioning

Abstract

While it has become axiomatic that mechanical signals promote in vitro engineered tissue formation, the underlying mechanisms remain largely unknown. Moreover, efforts to date to determine parameters for optimal extracellular matrix (ECM) development have been largely empirical. In the present work, we propose a two-pronged approach involving novel theoretical developments coupled with key experimental data to develop better mechanistic understanding of growth and development of dense connective tissue under mechanical stimuli. To describe cellular proliferation and ECM synthesis that occur at rates of days to weeks, we employ mixture theory to model the construct constituents as a nutrient-cell-ECM triphasic system, their transport, and their biochemical reactions. Dynamic conditioning protocols with frequencies around 1 Hz are described with multi-scale methods to couple the dissimilar time scales. Enhancement of nutrient transport due to pore fluid advection is upscaled into the growth model, and the spatially dependent ECM distribution describes the evolving poroelastic characteristics of the scaffold-engineered tissue construct. Simulation results compared favorably to the existing experimental data, and most importantly, distinguish between static and dynamic conditioning regimes. The theoretical framework for mechanically conditioned tissue engineering (TE) permits not only the formulation of novel and better-informed mechanistic hypothesis describing the phenomena underlying TE growth and development, but also the exploration/optimization of conditioning protocols in a rational manner.

Appendix

1.1 Non-dimensionalization of the governing equations

We non-dimensionalize the dependent variables of the governing equations of the growth model—for that, we employ 5 times the exterior oxygen concentration $5\rho _\mathrm{o}^{\mathrm{ext}}$, initial cell seeding $\rho _\mathrm{c}^{\mathrm{seed}}$, and maximum ECM density$\rho _\mathrm{m}^{\max }$ resulting in $\tilde{\rho }_\mathrm{o} ={\rho _\mathrm{o} }/{5\rho _\mathrm{o}^{\mathrm{ext}}}, \tilde{\rho }_\mathrm{c} ={\rho _\mathrm{c} }/{\rho _\mathrm{c}^{\mathrm{seed}}}$, and $\tilde{\rho }_\mathrm{m} ={\rho _\mathrm{m} }/{\rho _\mathrm{m}^{\max } }$, respectively. The specific form of the governing equations does not change, but upon non-dimensionalization, constants present in the model equations are updated accordingly (and their non-dimensionalized forms and values are recorded in Table 2). Non-dimensional forms of the governing Eqs.  (2.1), (2.3), and (2.5) will have tildes in all their members and are

$$\begin{aligned}&\dot{{\tilde{\rho }}}_\mathrm{o} +\nabla \cdot (-D_\mathrm{o} \nabla \tilde{\rho }_\mathrm{o} )+\hat{\mathbf{v}}\cdot \nabla \tilde{\rho }_\mathrm{o} =\tilde{\rho }_\mathrm{c} \frac{\tilde{Q}_\mathrm{o} \tilde{\rho }_\mathrm{o} }{\tilde{\rho }_\mathrm{o}^\mathrm{h} +\tilde{\rho }_\mathrm{o} } \end{aligned}$$

(6.1)

$$\begin{aligned}&\dot{{\tilde{\rho }}}_\mathrm{c} +\nabla \cdot \left[ {-D_\mathrm{c} \nabla \tilde{\rho }_\mathrm{c} +\tilde{\chi }\tilde{\rho }_\mathrm{c} \left( {\tilde{\rho }_\mathrm{c}^{\max } -\tilde{\rho }_\mathrm{c} } \right) \,\nabla \tilde{\rho }_\mathrm{o} } \right] \nonumber \\&\quad =\left[ {\tilde{k}^{\mathrm{p}}\frac{(\tilde{\rho }_\mathrm{o} )^{n}}{\tilde{\rho }_\mathrm{o}^\mathrm{c} +(\tilde{\rho }_\mathrm{o} )^{n}}\hbox {e}^{-\tilde{k}_\mathrm{i} \tilde{\rho }_\mathrm{o} }-\tilde{k}^{\mathrm{a}}} \right] \tilde{\rho }_\mathrm{c} \left( {\tilde{\rho }_\mathrm{c}^{\max } -\tilde{\rho }_\mathrm{c} } \right) \end{aligned}$$

(6.2)

$$\begin{aligned}&\dot{{\tilde{\rho }}}_\mathrm{m} =\tilde{k}_\mathrm{m} \tilde{\rho }_\mathrm{c} \left( {1-\tilde{\rho }_\mathrm{m} } \right) \nonumber \\&\qquad =\left\{ {{\begin{array}{ll} 0 , &{}\quad {\tilde{\rho }_\mathrm{o} (\mathbf{x},t)\le \tilde{\rho }_\mathrm{o}^s } \\ {\tilde{k}_\mathrm{m} (\tilde{\rho }_\mathrm{o} -\tilde{\rho }_\mathrm{o}^m )\tilde{\rho }_\mathrm{c} \left( {1-\tilde{\rho }_\mathrm{m} } \right) ,} &{}\quad \quad {\tilde{\rho }_\mathrm{o} (\mathbf{x},t)>\tilde{\rho }_\mathrm{o}^s } \\ \end{array} }} \right. \end{aligned}$$

(6.3)

and their parameters are adjusted to their non-dimensional form (listed in Table 2). As we will always refer to their non-dimensional form, tildes on the dependent variables and parameters will not be carried.

1.2 ECM-dependent diffusivity and oxygen-dependent cellular proliferation/apoptosis kinetics

We employ $\tilde{\rho }_\mathrm{m} ={\rho _\mathrm{m} }/{\rho _\mathrm{m}^{\max } }\in [0,1]$, a non-dimensional measure of ECM, to modulate oxygen diffusivity as the construct grows and evolves over time: $\tilde{\rho }_\mathrm{m} \rightarrow 0, D_\mathrm{o} (\tilde{\rho }_\mathrm{m} )\rightarrow D_\mathrm{o}^{\tilde{\rho }_\mathrm{m} =0}$ represents the initial diffusivity of oxygen in the scaffold devoid of ECM material, whereas $\tilde{\rho }_\mathrm{m} \rightarrow 1$, $D_\mathrm{o} (\tilde{\rho }_\mathrm{m} )\rightarrow D_\mathrm{o}^{\tilde{\rho }_\mathrm{m} =1}$ corresponds to the diffusivity of oxygen in native tissue or in engineered tissue at the limiting steady state of cultivation at a balance of production, degradation, and incorporation of collagen. Wood et al. (2002) have found that effective diffusivity in cellular systems is somewhat insensitive to the detailed geometric structure of the system, whereas it is primarily influenced by the volume fractions of the extracellular and intracellular spaces. Wood et al. (2002), Wood and Whitaker (2000) have shown that the analytical solution obtained from Chang’s unit cell (1983), composed of a periodic array of spheres and accounting for diffusion in the extracellular and intracellular spaces and finite mass transport across the cellular membrane, provides reasonable estimates of the experimentally measured effective diffusivity of several solutes in diverse biological tissues. We specify the ECM-dependent oxygen diffusivity as a function of $\tilde{\rho }_\mathrm{m}$:

$$\begin{aligned} \frac{D_\mathrm{o} (\tilde{\rho }_\mathrm{m})}{D_\mathrm{o}^{\tilde{\rho }_\mathrm{m} {=}0} }\,{=}\,\frac{3\kappa -2(1\,{-}\,\tilde{\rho }_\mathrm{m} )(\kappa {-}\,1)\,{+}\,2(1\,{-}\,\tilde{\rho }_\mathrm{m} )(3\tilde{\rho }_\mathrm{m} )^{-1}\kappa \psi }{3{+}\,(1\,{-}\,\tilde{\rho }_\mathrm{m} )(\kappa \,{-}\,1)\,{+}\,(2\,{+}\,\tilde{\rho }_\mathrm{m} )(3\tilde{\rho }_\mathrm{m} )^{-1}\kappa \psi }, \end{aligned}$$

(6.4)

where $\kappa$ and $\psi$ are non-dimensional material parameters, the former is a ratio of the diffusivities in extracellular and intracellular media and the latter weights transmembrane transport properties with extracellular diffusion. If membrane resistance is negligible, Eq. (6.4) recovers Maxwell’s solution of heat conduction in a heterogeneous two-phase system (Wood and Whitaker 2000). At $\tilde{\rho }_\mathrm{m} \rightarrow 1, D_\mathrm{o} (\tilde{\rho }_\mathrm{m} )\rightarrow D_\mathrm{o}^{\tilde{\rho }_\mathrm{m} =1}$ and Eq. (6.4) results in

$$\begin{aligned} \frac{D_\mathrm{o}^{\tilde{\rho }_\mathrm{m} =1}}{D_\mathrm{o}^{\tilde{\rho }_\mathrm{m} =0} }=\frac{3\kappa }{3+\kappa \psi } \end{aligned}$$

(6.5)

and from which specific values for parameters $\kappa$ and $\psi$ can be estimated.

The oxygen-dependent cellular proliferation/apoptosis rate $k_\mathrm{c} =k_\mathrm{c} (\rho _\mathrm{o} )$ describes the phenomena of growth inhibition and even apoptosis at very high oxygen toxic concentrations, and at very low oxygen anoxic concentrations, and is given by

$$\begin{aligned} k_\mathrm{c} =k_\mathrm{c} (\rho _\mathrm{o} )=k^{\mathrm{p}}\frac{(\rho _\mathrm{o} )^{n}}{\rho _\mathrm{o}^\mathrm{c} +(\rho _\mathrm{o} )^{n}}\mathrm{e}^{-k_\mathrm{i} \rho _\mathrm{o} }-k^{\mathrm{a}} \end{aligned}$$

(6.6)

The first term of Eq. (6.6) is based on Moser’s growth kinetics (driven by exponent n) and Aiba et al. inhibition kinetics (driven by the exponential contribution on the first term) to describe the inability for proliferate at low and at high oxygen concentrations, respectively. Kinetic constant $k^\mathrm{p}$ is a maximum proliferation rate, $\rho _\mathrm{o}^\mathrm{c}$ is the oxygen concentration for half-maximal proliferation rate, and exponent n is the anoxic inhibition coefficient for the Moser-type growth term. Kinetic constant $k_\mathrm{i}$ is the oxygen’s toxicity inhibition rate for Aiba-type inhibition term, and $k^\mathrm{a}$ the maximum apoptosis rate. Practically, Eq. (6.6) results in $k_\mathrm{c} <0$ for low and high concentrations of oxygen, and in $k_\mathrm{c} >0$ for a region in between where oxygen conditions allow proliferation to occur (Fig. 2).

1.3 Poroelasticity formulation

Let us consider a mixture composed of a porous solid saturated with fluid. As the mixture deforms, particles of the constituents that co-occupied a point move along different paths and this allows the description of the diffusion or flow of fluid with respect to the solid. We employ the u-p classical formulation of finite poroelasticity (Coussy and Coussy 2004). The dependent variables are the displacement field of the poroelastic mixture, denoted by u, and the pore pressure, denoted by $p^\mathrm{f}$ (Fig. 1). Let the elementary volume of the porous mixture in the reference and current configurations be represented by

$$\begin{aligned}&\hbox {d}V = \hbox {d}V_\mathrm{f} +\hbox {d}V_\mathrm{s} , \end{aligned}$$

(6.7)

$$\begin{aligned}&\hbox {d}v =\hbox {d}v_\mathrm{s} +\hbox {d}v_\mathrm{f} , \end{aligned}$$

(6.8)

which is made up of volume occupied by the solid constituent ($\hbox {d}V_\mathrm{s}$ and $\hbox {d}v_\mathrm{s}$) and the fluid saturated void space ($\hbox {d}V_\mathrm{f}$ and $\hbox {d}v_\mathrm{f}$). The porosity of the medium $\varphi$ is the ratio of the volume of fluid to the total volume of a representative volume element, i.e., $\varphi ={\hbox {d}v_\mathrm{f} }/{\hbox {d}v}$ and $\varphi _0 ={\hbox {d}V_\mathrm{f} }/{\hbox {d}V}$ in the current and reference configurations, respectively, and related by

$$\begin{aligned} \varphi =1-J_\mathrm{f} J^{-1}(1-\varphi _0 ), \end{aligned}$$

(6.9)

where $J_\mathrm{f} ={\hbox {d}v_\mathrm{f} }/{\hbox {d}V_\mathrm{f} }$ is the change in volume of void space, and $J={\hbox {d}v}/{\hbox {d}V}$ is the change in volume of the porous mixture. Balance of mass of fluid is given by

$$\begin{aligned} \frac{1}{J}\frac{\hbox {d}}{\hbox {d}t}(J\varphi \rho _\mathrm{f} )+\nabla \cdot (\rho _\mathrm{f} \varphi {\mathbf{v}_\mathrm{f}} )=0, \end{aligned}$$

(6.10)

where $\mathbf{v}_\mathrm{f}$ is the average velocity of the fluid relative to the solid (seepage velocity or Darcy’s velocity). The total stress acting at a point, T, is assumed to be composed by an average pressure stress in the fluid, $p^\mathrm{f}$ and the effective stress in the solid $\mathbf{T}^{s}$, and is defined by

$$\begin{aligned} \mathbf{T}=-p^\mathrm{f}\mathbf{1}+\mathbf{T}^\mathrm{s} \end{aligned}$$

(6.11)

Balance of linear momentum of the mixture in the absence of body forces and restricting the analysis to quasi-static motions is given by

$$\begin{aligned} \nabla \cdot \mathbf{T}=\mathbf{0} \end{aligned}$$

(6.12)

The constitutive behavior of pore fluid flow is governed by Darcy’s law. Darcy’s law states that the volumetric flow rate of the fluid through a unit area is proportional to the negative of the gradient of pore pressure, i.e.,

$$\begin{aligned} \varphi \mathbf{v}_\mathrm{f} =-\mathbf{K}\nabla ({p_\mathrm{f} }/{\rho _\mathrm{f} )} \end{aligned}$$

(6.13)

where K is the permeability of the poroelastic mixture (units of length/time). We consider that permeability is isotropic, i.e., $\mathbf{K}=k\mathbf{1}$, with k the scalar permeability. The effective stress in the solid is assumed to follow classical hyperelasticity. The Cauchy stress tensor in the solid $\mathbf{T}^\mathrm{s}$ relates to the first Piola–Kirchhoff stress tensor P by

$$\begin{aligned} \mathbf{T}^\mathrm{s}=J^{-1}\mathbf{PF}^{\mathrm{T}}, \end{aligned}$$

(6.14)

where F is the deformation gradient computed from the gradient of displacements as

$$\begin{aligned} \mathbf{F}=\mathbf{1}+\nabla \mathbf{u}, \end{aligned}$$

(6.15)

and $J={\hbox {d}v}/{\hbox {d}V}=\det \mathbf{F}$. Finally, P is derivable from a stored energy function $W=W(\mathbf{F})$ through $\mathbf{P}={\partial W}/{\partial \mathbf{F}}$.

1.4 Initial and boundary value problems

The rectangular constructs have dimensions $L=25\,\hbox {mm}, d=5\,\hbox {mm}$, and $h=1\,\hbox {mm}$, and are modeled as a two-dimensional domain $\Omega$ along their length under the state of plane strain—it is implicitly assumed that $\Omega$ represents the middle longitudinal cross section of the construct and that changes in longitudinal sections are not appreciable at least far from the boundary effects existing near their side edges. To aid the specification of the boundary conditions, we will consider a (x, y) Cartesian system with the axes x and y aligned with the length and the thickness of the construct, respectively. The origin is defined in the middle of the construct; thus, $\Omega =\{{-L}/2<x<L/2,{-h}/2<y<h/2\}$. The fixed post is the boundary at $\Gamma _F =\{x={-L}/2,{-h}/2<y<h/2\}$, whereas the moving post is represented by the boundary $\Gamma _M =\{x=L/2,{-h}/2<y<h/2\}$. Top and bottom surfaces are defined by $\Gamma _T =\{{-L}/2<x<L/2,y={-h}/2\}$ and $\Gamma _B =\{{-L}/2<x<{-L}/2,y=h/2\}$ (Fig. 3a).

Initial and boundary conditions must be specified to particularize the solution of governing equations of the TE growth model [Eqs. (2.1), (2.3), (2.5)]. Equation (2.5) governing matrix production does not necessitate boundary conditions, as it is a reaction-only equation. Oxygen concentration is imposed on top and bottom surfaces (Dirichlet) with the exterior concentration

$$\begin{aligned} \left. {\rho _\mathrm{o} (\mathbf{x},t)} \right| _{\mathbf{x}\in \Gamma _{T,} \Gamma _B } =\rho _\mathrm{o}^{\mathrm{ext}} , \end{aligned}$$

(6.16)

which is assumed to be constant along the length and not time varying. No-flux boundary conditions (Neumann) are imposed on the cellular phase

$$\begin{aligned} \left. {\nabla \rho _\mathrm{c} (\mathbf{x},t)} \right| _{\mathbf{x}\in \Gamma _{T,} \Gamma _B } =\mathbf{0}, \end{aligned}$$

(6.17)

with the meaning that cells will not leave the construct through its top and bottom boundaries. Lastly, the two ends of the construct are considered impervious for both nutrients and cells as these portions of the construct are connected to the rigid metal posts, and no flux is considered to occur across these boundaries. As for initial conditions, we consider that nutrient and cells are homogeneously distributed at initial time, with exterior oxygen concentration (i.e., no oxygen tension), initial cell seeding, and no ECM,

$$\begin{aligned}&\rho _\mathrm{o} (\mathbf{x},0)=\rho _\mathrm{o}^{\mathrm{ext}} , \nonumber \\&\rho _\mathrm{c} (\mathbf{x},0)=\rho _\mathrm{c}^{\mathrm{seed}} , \nonumber \\&\rho _\mathrm{m} (\mathbf{x},0)=0. \end{aligned}$$

(6.18)

Governing equations of poroelasticity (6.10) and (6.12) require boundary conditions to be solved. Similarly as before, we consider that the two ends of the construct are impervious to fluid flow,

$$\begin{aligned} \left. {\mathbf{v}_\mathrm{f} (\mathbf{x},t)} \right| _{\mathbf{x}\in \Gamma _F ,\Gamma _M} =\mathbf{0}, \end{aligned}$$

(6.19)

whereas top and bottom surfaces allow fluid movement in and out the poroelastic body. The exterior fluid pressure is imposed in the pores, i.e.,

$$\begin{aligned} \left. {p_\mathrm{f} (\mathbf{x},t)} \right| _{\mathbf{x}\in \Gamma _T ,\Gamma _B } =p^{\mathrm{ext}}. \end{aligned}$$

(6.20)

All edges of the construct are traction-free, and displacement is imposed at two corners of the construct to describe the effects of the fixed and moving posts, $\mathbf{x}_F =({-L}/2,{-h}/2)$ and $\mathbf{x}_M =(L/2,{-h}/2)$, respectively,

$$\begin{aligned}&\mathbf{u}(\mathbf{x}_F ,t)=\mathbf{0}, \nonumber \\&\mathbf{u}(\mathbf{x}_M ,t)=-u_M (t)\mathbf{e}_y , \end{aligned}$$

(6.21)

where $u_M (t)$ is a periodic function describing the movement of the post along the length direction. For computational ease, we consider the motion to be sufficiently smooth, with period T and amplitude $u_M^{\max }$ (Fig. 3b), and the following piecewise form was chosen

$$\begin{aligned} u_M (t)=\left\{ {{\begin{array}{l} 0, \quad {0\le t<\gamma } \\ {\frac{u_M^{\max } }{2}\left[ {1+\sin \left( {\frac{2\pi t}{T-4\gamma }-\frac{\pi }{2}\frac{T}{T-4\gamma }} \right) } \right] ,}\\ \quad {\gamma \le t<{-\gamma +T}/2} \\ {u_M^{\max } },\quad {{-\gamma +T}/2\le t<{\gamma +T}/2} \\ {\frac{u_M^{\max } }{2}\left[ {1+\sin \left( {\frac{2\pi t}{T-4\gamma }-\frac{\pi }{2}\frac{T+8\gamma }{T-4\gamma }} \right) } \right] ,}\\ \quad {{\gamma +T}/2\le t<-\gamma +T} \\ 0, \quad {-\gamma +T\le t<T} \\ \end{array}}} \right. \end{aligned}$$

(6.22)

where $\gamma \ll T$ is a short threshold time during which the initial and max deflected positions are held to ensure sufficient smoothness of the motion. The maximum deflection of the construct, if pure bending would be assumed (i.e., the deformed shape would lay in a circular arc of radius R and spanning angle $2\theta$), is given by $\delta _{\max } =R(1-\cos \theta )$ where R and $\theta$ are the solutions of the following system of equations: $2\theta R=L$ and $R\sin \theta =L-u_M^{\max }$. In order to achieve a maximum deflection of$\delta _{\max } =6.36\,\hbox {mm}$, the deflection employed by Engelmayr et al. (2005, 2006), the maximum displacement of the moving post is $u_M^{\max } =-5.00\,\hbox {mm}$.

We neglect inertial effects and restrict our analysis to quasi-static motions; therefore, initial conditions on mixture velocity $\dot{\mathbf{u}}(\mathbf{x},0)$ are not required. Pore pressure is considered to be homogeneously distributed and equal to the exterior pressure at initial times

$$\begin{aligned} p_\mathrm{f} (\mathbf{x},0)=p^{\mathrm{ext}}, \end{aligned}$$

(6.23)

which with Darcy’s law results in fluid velocity initially zero everywhere, i.e., $\mathbf{v}_\mathrm{f} (\mathbf{x},0)=0$.

1.5 Results of interest for experimental validation

Quantities of interest are the evolution of the distribution of cells and ECM over time, i.e., $\rho _\mathrm{c} (\mathbf{x},t)$ and $\rho _\mathrm{m} (\mathbf{x},t)$ . Whereas the control–case results in a one-dimensional problem (due to the particular form of the boundary conditions, the solution is independent of x), mechanical conditioning introduces a degree of inhomogeneity in the equation governing oxygen transport and consumption. Results will be presented either in the form of profiles at a given time t along the thickness of the scaffold at the mid-plane plane $x=0,\,y\in [{-h}/2,h/2]$ and contours over the entire length-section $\Omega$ (both corresponding to immunohistochemistry images of construct sections), and lastly, a single scalar corresponding to the total mass in the whole construct (corresponding to the assays conducted experimentally to the entire constructs) obtained with

$$\begin{aligned}&M_\mathrm{c} (t)=\int \limits _\Omega {\rho _\mathrm{c} (\mathbf{x},t)\hbox {d}v} \nonumber \\&M_\mathrm{m} (t)=\int \limits _\Omega {\rho _\mathrm{m} (\mathbf{x},t)\hbox {d}v} \end{aligned}$$

(6.24)

The temporal evolution of the flexural rigidity of the construct, which considering no stiffening effect due to deformation, is computed at the middle section with

$$\begin{aligned} EI(t)=\frac{E_\mathrm{s} h^{3}d}{12}+\bar{{E}}_\mathrm{m}^0 d\int _{-h/2}^{h/2} {\left. {y^{2}\rho _\mathrm{m} (\mathbf{x},t)} \right| } _{\mathbf{x}=(0,y,0)} \,\hbox {d}y \end{aligned}$$

(6.25)

Lastly, the temporal evolution of the force necessary to deform the constructs is obtainable from the poroelastic solutions.

1.6 Parameter estimation

Exterior oxygen concentration of 20 % $\hbox {O}_{2}$ is reported by Zhao et al. (2005) and Pathi et al. (2005), who have employed a similar boundary condition in their simulations, as $\rho _\mathrm{o}^{\mathrm{ext}} =2.10\times 10^{-7}\hbox {mol}\,\hbox {cm}^{-3}$. We consider experimental data obtained by Lis et al. (Lis et al. 1987) on the collagen content of the human mitral valve to set $\rho _\mathrm{m}^{\max } =108.3\,\hbox { mg}/\hbox {g ww}$—the authors have obtained a leaflet composition of 1.4 g wet weight, 81 % of the wet weight composed by water, and 57 % of dry weigh composed by collagen. Initial cell seeding in TE is usually conducted at sufficiently high values as similar to native tissue cell content, and during the initial periods of incubation, cell content decreases substantially. Engelmayr et al. (2006) reported that the DNA content measured for native juvenile ovine PV leaflets as $701\pm 42\, {\upmu } \hbox {g}/\hbox {g ww}$. Assuming DNA content of 7.6 pg/cell and a wet tissue density of $1\hbox { g}/\hbox {cm}^{3}$, cell density in native tissue is set as $\rho _\mathrm{c}^{\max } =92.2\times 10^{6}\,\hbox { cell}/\hbox {cm}^{3}$. On their experiment with VSMCs, the DNA content of $536\pm 38\,{\upmu }\,\hbox {g}/\hbox {g ww}$ after 30-h seeding period (Engelmayr et al. 2005)—thus, we set seeding density as $\rho _\mathrm{c}^{\mathrm{seed}} =70.5\times 10^{6}\,\hbox { cell}/\hbox {cm}^{3}$.

To characterize the Michaelis–Menten consumption kinetics of oxygen [cf. Eq. (2.2)], we employ constants used by Pathi et al. (2005) and set $Q_\mathrm{o} =1.25\times 10^{-17}\,\hbox { mol}_{\mathrm{O}_{2}} \,\hbox {cell}^{-1}\,\hbox {s}^{-1}$ and $\rho _\mathrm{o}^\mathrm{h} =5\,\% \hbox { O}_2 =1.105\times 10^{-8}\,\hbox { mol}_{{\mathrm{O}}_{2}}\, \hbox {cm}^{-3}$. Obradovic et al. (2000) considered that diffusivity of oxygen in TE scaffolds to be half of the diffusivity in water; Galban and Locke (1999a, b) considered that diffusivity in tissue to be 10 times lower than diffusivity in water—we set $D_\mathrm{o}^{\tilde{\rho }_\mathrm{m} =0} =1.5\times 10^{-5}\,\hbox { cm}^{2}\,\hbox {s}^{-1}$ and $D_\mathrm{o}^{\tilde{\rho }_\mathrm{m} =1} ={D_\mathrm{o}^{\tilde{\rho }_\mathrm{m} =0}}/{10}$. To fully characterize the dependence of diffusivity on existing ECM [cf. Eq. (6.4)], we will further assume that transmembrane transport is 10 times higher than intracellular diffusivity; thus, we obtain $\kappa =0.43$ and $\psi =23.07$. Relationship (6.4) is only mildly nonlinear (not shown), and a parametric study on its parameters has shown it to be somewhat insensitive to $\psi$ (the ratio of transmembrane transport to extracellular diffusivity) and mainly governed by $\kappa$ (the ratio of extracellular to intracellular diffusivities).

Cellular proliferation/apoptosis and extracellular matrix production rate constants are determined/estimated with the scant data existing in the literature. Galban and Locke (1999a, b) have conducted a parameter analysis on several growth parameters using a modified Contois growth kinetics function and have obtained cellular growth and death rates (with respect to dimensionless cellular volume fraction in between 0 and 1) on the orders of magnitude of $\pm 1\times 10^{-5}\,\hbox { s}^{-1}$; using a similar formalism, Pathi et al. (2005) and Zhao et al. (2005) estimated a maximum proliferation rates of $5.99\times 10^{-5}$ and $3.875\times 10^{-7}\,\hbox { s}^{-1}$, respectively—notwithstanding, it must be stressed that these proliferation and death rates are inherent to Galban and Locke’s framework, but can be employed to estimate the oxygen-dependent cellular proliferation/apoptosis rate function $k_\mathrm{c} (\rho _\mathrm{o})$ [Eq. (6.6); Fig. 3]. Experimental data and rate constants for the VSMCs employed by Engelmayr et al. (2005) do not exist; thus, we assume apoptosis rate $\tilde{k}^\mathrm{a}=1\times 10^{-5}\hbox { s}^{-1}$ and maximum proliferation rate $\tilde{k}_\mathrm{c}^{\max } =1.24\times 10^{-5}\,\hbox {s}^{-1}$, both specified with respect to changes in non-dimensional variable $\tilde{\rho }_\mathrm{c}$ and in agreement with orders of magnitude proposed by Galban and Locke (1999a, b). In order to obtain this behavior, we set constants in Eq.  (6.6) as: proliferation rate $\tilde{k}^\mathrm{p}=2.6\times 10^{9}$, inhibition kinetics rate $\tilde{k}^{i}=14$, Moser’s growth kinetics parameters $n=5$ and $\tilde{\rho }_\mathrm{o}^\mathrm{c} =0.45$ . Units of these parameters involve multiple powers of $\hbox {[mol}_{\mathrm{O}_{2}} /\hbox {cm}^{3}]$ because of exponent n, and are of much lesser relevance than the dimensionality of the full expression of $\tilde{k}_\mathrm{c} ,\, [\hbox {s}^{-1}]$. Engelmayr et al. (2005) have observed a marked region of non-ECM producing cells inside the construct, and we employed solutions of the oxygen governing equations with values of oxygen consumption kinetics and oxygen diffusivity above to estimate ECM production threshold of $\tilde{\rho }_\mathrm{o}^\mathrm{s} ={\rho _\mathrm{o}^s }/{(5\rho _\mathrm{o}^{\mathrm{ext}} })=16\,\% \hbox { O}_2$. Obradovic et al. (2000) have determined a GAG production rate of $2.3\,\% \hbox {ww}\,\hbox {day}^{-1}\,\hbox {mM}_{O_2 }^{-1} \times [10^{5}\hbox { cell}/\hbox {mm}^{3}]^{-1}$ occurring in cartilaginous TE experiments conducted by Vunjak-Novakovic and co-workers (Freed et al. 1998, 1994; Obradovic et al. 1999)—we employ such finding to estimate the ECM production rate $k_\mathrm{m}^0 =2.66\times 10^{-9}\hbox { mg}/\hbox {g ww}\,[\hbox {mol}_{\mathrm{O}_{2}} /\hbox {cm}^{3}]^{-1}\,[\hbox {cell}/\hbox {cm}^{3}]^{-1}\, \hbox {s}^{-1}$.

Engelmayr et al. (2005) reported Young’s modulus of the scaffolds as $E_\mathrm{s} =174\pm 22\hbox { kPa}$, with initial porosity of $\varphi _0^{\tilde{\rho }_\mathrm{m} =0} =0.96$, and that they are highly compressible—therefore, we set Poisson’s ratio of scaffold of ${\upnu }_\mathrm{s} =0.1$ and employ relationships $\mu =E/2(1+\nu )$ and $K=E/3(1-2\nu )$ to determine shear and bulk moduli of scaffold as $\mu _\mathrm{s} =39.54\hbox { kPa}$ and $K_\mathrm{s} =72.5\hbox { kPa}$, respectively, characterizing the initial condition of the two-parameter family of compressible neo-Hookean material describing the evolving poroelastic solid when no ECM is present. Engelmayr et al. (2005) have also observed a markedly linear relationship between collagen concentration and Young’s modulus of TE constructs (cf. Fig. 8 of Engelmayr et al. 2005)—although their range of observation is restricted to the initial stages of TE growth when a still small amount of ECM is present (up to a maximum of $\rho _\mathrm{m} \approx 1.5\hbox { mg/g ww}$ at 9 weeks of culture of TEHVs), it allows for parameter $\bar{{E}}_\mathrm{m}^0$ appearing in Eq. (2.14) to be obtained from its slope and set to $\bar{{E}}_\mathrm{m}^0 =0.8881\hbox { kPa}\,[{\upmu } \hbox {g}/\hbox {g ww}]^{-1}$. Finally, we assume that engineered tissue at its final steady-state stage of cultivation (when $\tilde{\rho }_\mathrm{m} \rightarrow 1$) is nearly incompressible, thus setting ${\upnu }_t =0.45$.

Variations of porosity and permeability with tissue evolution [set as linear relationships shown in Eqs. (2.7), (2.8)], require the specification of scaffold and tissue porosities and permeabilities and were estimated as follows: (i) Harrison and Massaro (1976) determined the water flux through a 2-mm-thick porcine aortic tissue due to a hydrostatic pressure gradient of 110 mmHg and have determined a hydraulic conductivity (corresponding to flux times thickness divided by pressure drop) of $7\times 10^{-13}\hbox { cm}^{4}\,\hbox {dyne}^{-1}\,\hbox {s}^{-1}$ from which we set permeability of tissue $k^{\tilde{\rho }_\mathrm{m} =1}=3.66\times 10^{-7}\hbox { cm}\,s^{-1}$; (ii) we employ the same reasoning as Galban and Locke (1999a, b) for the specification of the permeability of the scaffold (mainly due to lack of existing data) as 10 times greater than permeability of native tissue, i.e., $k^{\tilde{\rho }_\mathrm{m} =0}=10k^{\tilde{\rho }_\mathrm{m} =1}$; and finally, (iii) Lovich and Edelman (1996) have determined experimentally the fractional space in arterial media, adventitia, and myocardium, and we employ their result on arterial media to set porosity in tissue as $\varphi _0^{\tilde{\rho }_\mathrm{m} =1} =0.61$. The remaining parameter employed for the poroelastic formulation are fluid density, set to density of water $\rho _\mathrm{f} =1\hbox { g}/\hbox {cm}^{3}$, and external pressure, set to atmospheric pressure $p^{\mathrm{ext}}=101.325\hbox { kPa}$.

Lastly, under static conditions nutrient transport occurs due to diffusion only, i.e., $\hat{\mathbf{v}}=\mathbf{0}$ in Eq.  (2.1), and the framework, together with the parameter set above, is able to describe with sufficient accuracy the existing experimental data set (cross-sectional profiles, bulk measures, and bulk mechanical properties). When subjected to mechanical conditioning, the mechanism of nutrient convection plays a significant role in changing the behavior of the system. Time-lumping constant C present in Eq. (2.15) was set to $C=70$ such that the departure from the static solution represents the observed differences in between both regimes by Engelmayr et al. (2005) with the same parameter set. Although we have included in our model the effects of chemotaxis and of deformation-dependent ECM production (both in quantity and in quality), we have chosen to disregard these effects mainly due to the lack of critical experimental data to determine their impact. Consequently, chemotactic flux in Eq. (2.3) is not accounted, and matrix production rate and stiffness [Eqs. (2.6), (2.10)] are indifferent to deformation, i.e., $\bar{{k}}_\mathrm{m} (\mathbf{C})=0$ and $\bar{{w}}_\mathrm{m} (\mathbf{C}^{t}(s))=0$. Notwithstanding, it must be stressed that in general these effects might be relevant and significant in engineered tissue growth and that motivates our reasoning to explicitly account deformation-dependent ECM production rate and stiffness at the model development stage. The determination of these effects must be carried on with different experimental programs than the one conducted by Engelmayr et al. (2003, 2005, 2006) and are subject of our future studies.