A poroelastic mixture model of mechanobiological processes in biomass growth: theory and application to tissue engineering.

In this article we propose a novel mathematical description of biomass growth that combines poroelastic theory of mixtures and cellular population models. The formulation, potentially applicable to general mechanobiological processes, is here used to study the engineered cultivation in bioreactors of articular chondrocytes, a process of Regenerative Medicine characterized by a complex interaction among spatial scales (from nanometers to centimeters), temporal scales (from seconds to weeks) and biophysical phenomena (fluid-controlled nutrient transport, delivery and consumption; mechanical deformation of a multiphase porous medium). The principal contribution of this research is the inclusion of the concept of cellular "force isotropy" as one of the main factors influencing cellular activity. In this description, the induced cytoskeletal tensional states trigger signalling transduction cascades regulating functional cell behavior. This mechanims is modeled by a parameter which estimates the influence of local force isotropy by the norm of the deviatoric part of the total stress tensor. According to the value of the estimator, isotropic mechanical conditions are assumed to be the promoting factor of extracellular matrix production whereas anisotropic conditions are assumed to promote cell proliferation. The resulting mathematical formulation is a coupled system of nonlinear partial differential equations comprising: conservation laws for mass and linear momentum of the growing biomass; advection-diffusion-reaction laws for nutrient (oxygen) transport, delivery and consumption; and kinetic laws for cellular population dynamics. To develop a reliable computational tool for the simulation of the engineered tissue growth process the nonlinear differential problem is numerically solved by: (1) temporal semidiscretization; (2) linearization via a fixed-point map; and (3) finite element spatial approximation. The biophysical accuracy of the mechanobiological model is assessed in the analysis of a simplified 1D geometrical setting. Simulation results show that: (1) isotropic/anisotropic conditions are strongly influenced by both maximum cell specific growth rate and mechanical boundary conditions enforced at the interface between the biomass construct and the interstitial fluid; (2) experimentally measured features of cultivated articular chondrocytes, such as the early proliferation phase and the delayed extracellular matrix production, are well described by the computed spatial and temporal evolutions of cellular populations.

Introduction

It is nowadays a founding concept in cellular and molecular biology that cells are able to sense mechanical stimuli in their surrounding environment and produce a coordinate response. Such a process, defined as mechanotransduction (see, e.g., the recent review in [50]), plays important roles in several physiological processes such as cell motility, angiogenesis, bone formation and wound healing [76]. In this work, we present a mathematical approach for describing mechanotransduction processes involved in tissue growth. The proposed description, albeit very general, is applied to the scenario represented by tissue engineering. In this context, a better knowledge of the role of biomechanical cues can help in orchestrating a more effective artificial tissue growth. More in detail, our work is motivated by a specific tissue engineering application, artificial regeneration of articular cartilage. Briefly, cartilage cells (articular chondrocyte cells, ACCs) or other progenitor cells are seeded into polymeric scaffolds, possibly perfused by an interstitial fluid to force nutrient delivery. Cells are expected to duplicate and, above all, produce an increasing mass of ECM, forming cartilagineous neo-tissue. Cartilage tissue growth in engineered constructs had been already studied in a series of papers by Klisch and coauthors [30-32]. In this work, we enrich the description of the biophysical phenomena by introducing the conceptual framework developed in [43, 45, 46]. In these works, the isotropic/anistropic state of the cytoskeletal tension is shown to be responsible for triggering signalling transduction cascades which regulate functional cell behaviors related to proliferation and/or ECM secretion. Under the assumption of an isotropic strain-stress response, a uniform distribution of stress over the cell surface - stress due to the traction forces exerted by the cell on the surrounding environment - generates an “isotropic cytoskeletal tension state” in which the cell nucleus tends to maintain a roundish morphology (see Fig. 1a). Conversely, in an “anisotropic cytoskeletal tension state” the cell nucleus tends to elongate (see Fig. 1b). Macroscopically speaking, when the nucleus maintains a roundish morphology, ECM secretion is favoured, whereas an elongated nucleus favors cell duplication by dividision along a polarization axis represented by its longer axis itself. According to a finer biomolecular view, this process can be interpreted as due to the fact that the shape of the nucleus is known to regulate the porosity of its membrane and, through this, the import flow of specific transcription factors which regulate cell behavior (we refer to [43] and references therein for a more biologically detailed discussion of these complex processes).

Fig. 1

Concept of isotropic/anisotropic stress state of a cell: a a uniform distribution of traction forces (mediated by cell membrane integrin/cadherins) over the cell surface generates an “isotropic cytoskeletal tension state” in which the cell nucleus tends to maintain a roundish morphology, favoring ECM secretion; b an “anisotropic cytoskeletal tension state” the cell nucleus tends to elongate, favoring cell mitosis

In order to work in a continuum mechanics framework, we propose an extension of the above concepts to aggregates of cells. Figure 2 schematically illustrates how the mechanobiological conditions may affect and drive the fate of a colony of ACCs seeded in a 3D porous scaffold. When cells are first seeded in the scaffold, they form a thin layer covering the surface. Since the characteristic dimension of the local scaffold curvature is much larger than cell size, cells find themselves in a local planar condition (see Fig. 2a). Exerting adhesion forces on the scaffold surface, cells tend to assume a spread elongated shape, orienting themselves along a preferred polarization axis. According to the concept of “force isotropy”, this represents a condition which enhances the probability that the single cell enters into a proliferative state. This situation persists until all the pore surface is covered with cells (Fig. 2b). From this moment on, cells start to occupy the empty space of the pore (Fig. 2c, d). Cells in contact with other cells sense an isotropic stress state condition, which drives the cell towards a mature differentiated phenotype, characterized by ECM secretion (Fig. 2e).

Fig. 2

Various phases of tissue growth inside a scaffold pore: a seeding phase and cell polarization; b proliferation and formation of a monolayer; c, d formation of new construct layers; e ECM secretion

To translate the mechanobiological description of the processes illustrated in Fig. 2 into a mathematical model, we combine the poroelastic theory of mixtures and cellular population models. Using the model, we predict the spatial and temporal distribution of a biomass aggregate of ACCs and ECM under the biophysical assumption that stress state and oxygen (nutrient) tension act as main determinants of engineered culture evolution. The poroelastic theory of mixtures has already been proposed elsewhere to describe mechanobiological processes in growing tissues, possibly combined with a multiscale approach. In these formulations, tissue growth is represented as mass exchange between phases in a globally mass-conserving framework. Moreover, the assumption of linear strains is made, justified by the fact that the microscopic representative volume in which the volume averaging is performed does not evolve in time whereas the individual phases do. We refer for these approaches to the comprehensive discussion in [5], where several different techniques (effective medium theory, mixture theory, volume averaging and asymptotic (two-scale) homogenization) available to describe a poroelastic growing system are discussed and compared. Moreover, we refer to [54] for an example of the use of asymptotic homogenization techniques to develop a model for growing poroelastic media. As for cellular population models, they are used in several literature papers to describe the evolution of mixture components. We refer in particular to the works of [51] and [67], where multiple cellular populations are studied describing the exchange from one population to the other via a phenomenological representation. From this perspective, the principal contribution of our mechanobiological model is the inclusion of the concept of cellular “force isotropy” as a determinant of the passage from one pool of cells to the other (proliferating, ECM secreting o quiescent cells). A phenomenological indicator of the stress/strain state of the continuum construct is proposed, based on the on the norm of the deviatoric part of the total stress tensor computed via the poroelastic theory. Similarly to other models in tissue engineering applications (see, e.g., [14, 47], we also include the effect of nutrient (oxygen) availability by solving for it a transport-diffusion-reaction equations. Nutrient levels are supposed to be as well driving mechanisms in the cellular pool exchange. We use the model on a preliminary simplified geometrical one-dimensional setting to study the influence of the different parameters on the evolution of the construct composition. The extensive numerical simulations carried out under different working conditions show that fundamental roles are played by the maximum cell specific growth rate and by the mechanical boundary conditions at the interface between biomass construct and interstitial fluid.

The paper is organized as follows: in Sect. 2, we present the assumptions leading to the description of the biomass as a mixture; in Sect. 3, we discuss the kinematic laws for the growing biomass; in Sect. 4 we formalize the balance laws for the biomass; in Sects. 5 and 6 we discuss the proposed exchange pathways interconnecting biomechanical cues and cell population evolution along with the corresponding model; in Sect. 7, we present the 3D mathematical model with the boundary conditions, while in Sect. 8 we introduce the reduced 1D setting with the proposed stress indicator; in Sect. 9, we discuss the numerical approximation of the 1D model and in Sects. 10 and 11 we present the results of the numerical simulations and we carry out a comprehensive discussion. Eventually, in Sect. 12 we draw the conclusions and we present perspectives for future work.

Multiphase modeling

In this section we develop a mathematical model based on the representation of the ensemble of the growing cartilaginous biomass by the mixture theory. In this framework, equations are postulated for the balance of mass and momentum for each constituent and then for the entire mixture according to the following ideas:

the growing biomass is treated as a mixture composed by a multiphase solid mass and an interstitial fluid, the latter representing a fraction of the order of 65–80% in mass of the total biomass. The multiphase solid consists of ACCs and of ECM. ACCs are pooled in different populations according to their life cycle status (proliferative, ECM secreting, quiescent), as in the works of Sengers [66] and Ducrot [25];

the poroelasticity theory is used to model the interaction of deformation and fluid flow in the fluid-saturated porous, elastic solid [6, 22];

the kinematics of the solid phase of the mixture is based on an infinitesimal–deformation approach, including the effect on the stress field of biological growth, according to the formulation proposed by Klisch and co–authors [31-33];

the mass conservation balance for each single constituent and for the mixture are written according to the formulation introduced by Lemon and co-authors [35, 36] and extensively analyzed in [51, 70];

the mass exchange terms, including the rate of switch of cells from a population to the other, are tuned according to the nutrient level, the latter being itself an unknown of the problem, and to the stress state locally experienced by the mixture, which may drive cells into a certain functional behavior pool.

In the following, we use the term “phase” when we refer to the solid or to the fluid part of the mixture, while the term “component” is used to refer to any of the constituents of the solid phase (cell populations and ECM). When it is not necessary to distinguish between phase and component, we simply use the term “species”. The meaning of the subscripts used throughout the article is as follows: s=solid phase, fl=fluid phase, cells=cell component of the solid phase, ECM= extracellular component of the solid phase.

We let and t denote the space and time variables, respectively. We use the convention that the dependence of all variables and model parameters on and t is left understood except otherwise stated.

The geometrical configuration of the mixture is identified by the open bounded set ( unless otherwise specified). The domain does not evolve in time, rather, it is the amount of each species at a point that changes with t due to cell proliferation and matrix deposition. This is the precise sense of the concept of “growing mixture”. From now on, we denote by the space-time cylinder in which the TE problem is studied, being the final time of culture process.

Referring to Fig. 3, for all we associate with a generic point a fixed representative elementary volume (REV) (see [73]) and denote by its d-dimensional volume. Then, for each component of the growing mixture, we define the volume fractionas the time evolving ratio of the volume occupied by the i-th component in the REV and the volume of the REV itself. We also let

Fig. 3

A schematic view of the computational domain with a detailed view of a typical REV where the various phases and components of the growing mixture are identified

According to the biochemical hypothesis a), we consider three ACC populations: proliferating, ECM secreting state and quiescent state, denoted by the letters , and , respectively. Then, we haveEventually, we denote by the vector-valued function comprising the volume fractions of the fluid phase, the three cellular populations and the ECM.

The following assumptions on the mixture are also considered.

Assumption 1

(Fully saturated mixture) The mixture is fully saturated, i.e.

Relation (1c) is referred to as saturation condition [3, 51, 65] and excludes the possibility of the formation of voids or air bubbles inside the growing mixture.

Assumption 2

(Intrinsic incompressibility) All species constituting the growing mixture have the same (constant) mass density of the physiological interstitial fluid assimilated to water [3, 26, 30, 36, 51]. This is not, in general, equivalent to assuming that the whole mixture is incompressible (see [51] p. 629).

Assumption 3

(Closed mixture) The mixture is closed, this meaning that the system does not exchange mass with external mass sources or sinks [51].

Kinematics of the growing mixture

From now on, we denote “solid matrix” the collection of solid phase constituents, that is, cells and ECM. Then, we apply to the solid matrix the so–called intermingled mixture constraint [2], stating that all the solid matrix constituents experience the same overall motion. This hypothesis amounts to assuming the displacement and velocity vectors of each constituent to coincide with those of the solid matrix. Then, we denote by and the displacement and velocity at the time level t of any point of the solid component of the biomass, and by the associated infinitesimal deformation of the biomass volume surrounding the point at time t. The intermingled mixture constraint yields also the following relation The growth process of each mixture component (cellular growth and ECM secretion) is taken into account by introducing the following decomposition [33]where is the infinitesimal growth tensor associated with each solid constituent of the biomass, accounting for the amount and the spatial orientation of the newly deposited mass, and is the elastic accommodation tensor necessary to reinforce at each time level the continuity of the whole solid upon growth. Finally, we denote by the relative velocity [51, 65] of the fluid phase with respect to the solid phase in the biomass, being the velocity of the interstitial fluid. For notational brevity, from now on, we simply write and instead of and , respectively.

Balance laws for the deformable growing biomass

In this section we illustrate the set of conservation laws that model the mechanobiological processes regulating biomass tissue growth. For further discussion and analysis of mixture theory applied to tissue growth modeling we refer the reader to [11, 36, 51].

Mass balance for mixture components

The mass balance equation for the growing mixture is given by the following coupled system of PDEs in conservation form to be solved in : where is the flux matrix and is the net mass production rate for which the following constraint holds, due to Assumption 3 Equation (3b) represents a phenomenological description of the flux density of each species under the effect of convective transport due to the fluid and solid velocity, respectively. It is worth noting that in cartilage tissue growth, cells do not typically exhibit a significant diffusive motion, rather, they need a solid support for surviving and for developing their functional activities (property of anchorage-dependence [47]). For this reason in the present work we neglect the contribution to the flux density due the the diffusive transport, unlike in other applications where this term plays a significant role [37, 41, 49].

Momentum balance for mixture components

Under the assumption of negligible inertial terms and absence of body forces and volumetric fluid mass sources and sinks, the linear momentum balance equation for the solid and fluid phases of the growing mixture is expressed by the following PDEs in conservation form to be solved in : where is the effective stress tensor of the component of the solid phase of the mixture, is the pressure exerted by the fluid phase and is the identity tensor. The isotropic stress accounts for the coupling, typical of poroelasticity, between the flow of the fluid and the deformation of the solid matrix, and in particular describes the contribution to the stress due to the fluid pressure within the structure.

The quantities , , are the total stress tensors of the solid and fluid phases, while are the interphase forces [36]. As usual, we neglect the effective stress tensor of the fluid, meaning that we assume that the internal fluid viscosity is negligible compared with the friction between the fluid and the solid matrix [4, 26, 51]. For the mathematical characterization of the forces we refer to [36] and [51]. We observe that, for all and at all , it holds

Total mass and momentum balance for the growing biomass

The mass balance equation for the whole growing mixture is obtained by summing each component in system (3a), using (2c) and Assumptions 1 and 3: where is the composite velocity of the mixture (cf. Eq. (2.4) of [51]). A simple manipulation allows us to write Eq. (5b) asIn a similar manner, summing Eqs. (4) and using (4e), we get the total momentum equation where is the total stress in the mixture.

Mass balance for nutrient concentration

The mass balance system (3) for the solid and fluid phases of the growing mixture is accompanied by a corresponding continuity equation for the nutrient concentration (oxygen) that is transported throughout the growing on mixture by the interstitial fluid. This continuity equation is expressed by the following PDE in conservation form to be solved in : the interstitial fluid velocity being computed using (2c) asThe mathematical description of the oxygen diffusion coefficient adopted in this article is the so-called Maxwell model [74], that allows to account, in a volume-averaged sense, for the microscopic composition of the biomass. More precisely, we introduce the effective diffusion coefficientwhere and represent the nutrient diffusivity in the fluid and solid phase, respectively, while is the coefficient regulating local mass equilibrium between nutrient concentration in the solid and fluid phases (see [74] for a detailed discussion).

The time rate of oxygen consumed by the cellular populations is modeled by a generalized form of the Michaelis-Menten kinetics where , , is the nutrient consumption rate of the cellular population and is the half saturation constant. We refer to [63] and the literature cited therein for a similar treatment of the oxygen consumption term in the framework of a multi-phase growing mixture.

Mass exchange pathways

The production terms , cells,ECM, introduced in Eq. (3c) mathematically describe the mechanisms of addition and/or removal of mass for each species constituting the biomass growing mixture.

The exchange between the different functional cellular pools are supposed to be mediated by local population concentration, local stress state, local nutrient concentration and by natural decay times (see Fig. 4). To quantify the stress-mediated effect we proceed as follows. Let H(z) be the Heaviside function such that for and for , z being a real number. Then, the stress-state dependent effect is represented by , r and being an indicator of the isotropy or anisotropy of the local stress state and a threshold value, respectively. If , the local state of stress is considered as isotropic, in the other case is considered anisotropic. According to our mechanobiological picture, an anisotropic stress state enhances transition towards the proliferative state, whereas an isotropic stress state enhances transition towards the ECM secreting state. We use a similar approach to quantify the concentration-mediated effect by introducing the indicator , being a threshold concentration for cell activity. For notational brevity, we let and .

Fig. 4

Conceptual scheme of exchange pathways among cellular populations (generalization of Fig. 5.3 in [67])

The following exchange/production rate terms are considered:According to the exchange laws illustrated above, the production terms associated with cell populations are defined as:

mitotic proliferation: we suppose cells in pool n proliferate at a rate given by Relation (7a) is a phenomenological law in which the term keeps into account contact inhibition effects, while the term is a nutrient–dependent modulation (Monod law [21]), being the half-saturation constant and the maximum growth rate, respectively;

ECM production rate: we consider here glycosaminoglycan (GAG) as the main marker for ECM accumulation [47] and we assume a simple proportionality law between the total amount of ECM and the secretion rate of GAG. Following [47] we assume the ECM production rate to be given by In the above relation, is the volume of a single cell, the constant of proportionality accounts for the heterogeneous composition of cartilagineous ECM (water for 70–80% of its wet weight, collagen fibrils for 10–15% and GAG for 5%) [14], c is oxygen concentration and a growth factor. The last term model the fact that ECM synthesis attains its maximum value when no extracellular matrix is present because more space is available for matrix production. Then, as soon as sythesized matrix accumulates at each point of the biomass construct, the available space diminuishes until reaches a maximum value and matrix secretion ceases. Notice that in the description of GAG secretion, we are assuming that at the initial time level, biomass is constituted by a uniform layer of cells and matrix (see [57] for a similar approach). This corresponds to neglecting the very initial phase where the seeded cells proliferate and “pave” the scaffold wall, and is consistent with the mathematical fact that a continuum-based approach does not enable to reproduce the subcellular mechanisms that regulate the early mitotic process. These latter processes should be more properly described by treating seeded cells as individual units that behave according to cellular automata schemes [16, 18, 27, 28].

decay pathways: all cellular compartments may evolve into quiescent (absence of cell activity due to an insufficient oxygen intake [20, 67]) or apoptotic phases (cellular death). Quiescence occurs if nutrient concentration c falls below the critical level , whereas apoptotic phase is related to age dependent cell death [60]. The time rates of change between state () and the inactive states (quiescence or apoptosis) are and , respectively;

exchange rate between pools : a first contribution in the direction is regulated by the mitotic characteristic inverse time constant and take the form A second contribution is regulated by the probability rate that a cell in pool enters into pool , enhanced by the mechanical factor , giving the rate term

exchange rate between pools : the probability rate that a cell in pool enters into pool and the opposite for are mediated by the mechanical terms and , respectively, to signify that anisotropy favors proliferation while isotropy ECM secretion.

To conclude the mathematical description of mass exchange terms, we define the extracellular fluid production in such a way to satisfy Assumption 3 and, consistently, relation (3d)From a biophysical point of view this is equivalent to assuming that mass exchanges occur only among cells/ECM and fluid, meaning that dead cells and degrading ECM are deteriorated into extracellular fluid, and conversely that the latter is “consumed” whenever cells duplicate or secrete ECM [70]. From a computational point of view, relation (8e) allows us to eliminate the dependent variable and the corresponding mass balance equation from system (3a) as done in [70], Sect. 2.2, in such a way that the fluid volume fraction can be computed by simple post-processing as

Bio-mechanical models for the deformable growing biomass

In this section we provide a mathematical description of the mechanobiological phenomena involving growth processes (cell duplication and ECM secretion). To this purpose, we introduce suitable bio-mechanical models for the growth tensors in the decomposition (2b) by extending the theory developed in [32] and [33].

Growth laws

We propose the following definitions of the growth tensors: where the symbol represents the tensor dyadic product and , are growth coefficients for which a model equation is provided below. Eqns. (9) state that the mass increment of each mixture solid constituent is isotropically deposited for the , and components, while is accumulated along a specific polarization direction, identified by the unit vector , for proliferating cells.

Biophysical motivations support our choice of the growth laws (9). Firstly, according to the concept of “force isotropy” on the cell introduced in [45, 46], cells that occupy the bio-synthesizing compartment ( compartment) experience an isotropic adherence condition and consequently tend to assume a spherical shape (see Fig. 5, left) whereas cells that live in the proliferating compartment ( compartment) are subjected to an anisotropic adhesion state and tend to elongate (see Fig. 5, right). Secondly, according to the infinitesimal deformation growth theory developed in [32], the deformation of an infinitesimal sphere of biomass growing into an ellipsoid can be reasonably described by an anisotropic growth tensor, while the deformation of an infinitesimal sphere growing into a larger sphere can be characterized by a isotropic growth tensor. For the sake of simplicity, the infinitesimal growth tensor for the species and are supposed to be isotropic. The definition of in Eq. (9b) and the law for its time evolution is a delicate issue. In [8], is characterized according to the dynamics of the evolution of the cell cytoskeleton, which reorganizes itself according to its mechano–sensing mechanisms. A simplified version of the model proposed in [8], and also adopted in the present work, is represented by the choicewhere is the normalized eigenvector of the infinitesimal strain tensor , associated with the eigenvalue of largest module, which physically corresponds to the maximum principal dilatation of the biomass around such a point [68].

Fig. 5

Cellular level (top): pictorial representation of the isotropic/anisotropic adherence condition. Continuum level (bottom): isotropic/anisotropic biomass growth

The coefficients , , give a measure of the amount of mass of the cellular population of type deposited at time t at point . To determine these quantities, we proceed as in [32] and [33] and require the following growth continuity initial value problem to be satisfied for all and for :

The quantity represents the amount of mass of the cellular population deposited at time t at point per unit time and per unit reference mass. According to the general indications illustrated in Sect. 2.2.4 of [33], the growth laws are phenomenological equations that indirectly describe chemical processes responsible for growth and can be typically expressed as “synthesis” rate minus a “degradation” rate, that may include a mass conversion rate from one constituent of the mixture to another. Also, the constants that appear in a specific growth law may depend parameterically on biological factors such as, for example, the level of a specific growth factor. Thus, based on the description carried out in Sect. 5, we set , in such a way that the initial value problems that furnish the characterization of the growth coefficients become, for : and having used the identities and .

Constitutive equations for the mechanical and fluidsubsystems

We assume that cells and ECM behave like linear elastic solids, so that the effective stress tensors associated with the solid components of the biomass are defined as where and are the Lamé parameters of each component of the solid phase, , and are the growth strain tensors introduced in (9). More sophisticated constitutive models might be adopted [3, 24, 42, 51], but their use is beyond the scope of the present work which is mainly devoted to proposing a computationally feasible mechanobiological model of in vitro cartilage tissue growth. We assume the relative velocity in Eq. (5c) to be expressed by the Darcy law (see, e.g., [11] and references cited therein)where the isotropic permeability tensor is defined as in [65], being a friction coefficient. To provide a physically consistent characterization of we apply the classic Stokes theory for viscous drag to the biomass mixture and obtain and being cell radius and interstitial fluid dynamic viscosity, respectively, from which we get

The mathematical model in 3D

In this section we summarize the mathematical model for the multiphase mixture constituting the growing tissue. We refer to Table 1 for the definition of all model parameters and their quantitative value used in numerical simulations. In what follows, we denote by an open bounded set of representing the computational domain in which the biomass growth process physically takes place, and by the boundary of on which an outward unit normal vector is defined (see Fig. 6). In view of the definition of the boundary conditions to be supplied to the mechanobiological model, it is convenient to indicate by the set of the dependent variables of the problem. Then, for each we assume that is divided into two disjoint portions, denoted by and , such that , and where Dirichlet and Neumann boundary conditions are applied, respectively. By doing so, if the boundary conditions are for the model block of Sect. 4.3, if the boundary conditions are for the model block of Sect. 4.4 whereas if the boundary conditions are for the model block of Sect. 4.1. We notice that, in general, there is no geometrical relation between the various decompositions of the domain boundary belonging to the same typology of boundary condition. The sole requirement is that, for each considered dependent variable, the union of the Dirichlet and Neumann partitions is the whole domain boundary . For a similar treatment of the splitting of the domain boundary we refer to [26] and [7].

Table 1

Numerical values of model parameters used in the simulation tests

Symbol	Definition	Value	Units	References
c_0	O_2 concentration for t=0	5×10^-6	gcm^-3	This work
c_sat	O_2 saturation concentration	6.4×10^-6	gcm^-3	[14]
c_thr	O_2 threshold concentration	1.6×10^-6	gcm^-3	[38]
c_apo	O_2 apoptosis concentration	3.2×10^-7	gcm^-3	[38]
K_eq	O_2 local mass equilibrium coefficient	0.1	-	[13]
D_c,s	O_2 diffusivity in the solid phase	0.75×10^-5	cm^2s^-1	[13]
D_c,fl	O_2 diffusivity in the fluid phase	1×10^-5	cm^2s^-1	[13]
V_b	Inlet velocity of perfusion fluid	50×10^-4	cms^-1	[14]
T_b	Stress due to perfusion fluid	100	mPa	[14]
μ_fl	Fluid dynamic viscosity at 20^∘ C	1.002·10^-2	gcm^-1s^-1	[29]
R_n=R_v	O_2 consumption rate for n/v-cells	3.9×10^-8	g(cm^3s)^-1	[63]
R_q	O_2 consumption rate for q-cells	10^-8	g(cm^3s)^-1	This work
K_1/2	O_2 half saturation constant	3.2×10^-6	gcm^-3	[63]
β_A→B	Transition rate from state A to state B	10^-5	s^-1	This work
k_apo	Apoptosis transition rate	3.858×10^-7	s^-1	[63]
k_qui	Quiescence transition rate	3.858×10^-7	s^-1	This work
k_deg	ECM degradation rate	7.7×10^-7	s^-1	[71]
k_g0	Maximum specific cell growth rate	5.8×10^-6	s^-1	[63]
k_g1	“low” specific cell growth rate	1×10^-7	s^-1	This work
k_g2	“high” specific cell growth rate	1×10^-5	s^-1	This work
E	Expansion coefficient	20		[14]
k_GAG	GAG synthesis rate	8.61×10^-11	cm^6(cellsg)^-1	[14]
K_sat	Monod saturation constant	1.927×10^-6	gcm^-3	[15]
D_η	Cells and ECM diffusion coefficient	1×10^-9	cm^2s^-1	This work
λ_η	Cells and ECM Lamé’s parameter	5.1937×10^3	dynecm^-3	[64]
μ_η	Cells and ECM Lamé’s parameter	1.8248×10^3	dynecm^-3	[64]
ϕ_ECM,max	Maximum ECM volume fraction	0.1		This work
R_cell	Cell radius	5×10^-4	cm	This work
V_cell	Cell volume	5.236×10^-10	cm^3	This work
τ_m	Mitotic characteristic time	172800	s	[65]
K_ref	Reference permeability	1.67·10^-5	cm^3sg^-1	This work

Fig. 6

Schematic representation of the 3D computational domain. The boundary conditions for the mechanical block of the model are considered. In this example, whereas is nonvanishining only on a subset of

Having specified the geometrical setting of the biomass growth process, the mathematical model proposed in the present article consists of the following PDE subsystems:

Poroelastic IV-BVP for biomass: for given , find the solid displacement , the fluid friction velocity and pressure p such that the following equations are satisfied in : supplied with the following initial and boundary conditions: where is the initial mixture displacement in the domain whereas and are the given displacement and traction fields defined on the decomposition of the computational domain boundary in the mechanical block of the model.

Mass balance IV-BVP for nutrient concentration: for given , and , find the nutrient concentration c and the nutrient flux density such that the following equations are satisfied in : supplied with the following initial and boundary conditions: where is the initial nutrient concentration in the domain whereas and are the given concentration and diffusive nutrient flux defined on the decomposition of the computational domain boundary in the nutrient block of the model.

Mass conservation IV-BVP for cellular populations: for given , , and , find the volume fractions such that the following equations are satisfied in : where , with supplied with the following initial and boundary conditions: where is the initial distribution of cellular volume fraction in the domain. Condition (13j) amounts to assuming that no cellular flux is exchanged with the external environment during the biomass growth process.

The mechanobiological model in 1D

In this section we formulate the proposed mechanobiological model in a one-dimensional (1D) geometrical configuration. This is a first step toward the simulation of a realistic structure such as the 3D scaffolded bioreactor used in the experimental analysis discussed in [34]. The 1D-formulation is constructed by describing the biomass as a nonhomogeneous bar (fixed at one endpoint) subject to a uniaxial state of mechanical stress in such a way that each point of the bar undergoes the same deformation. Then we consider the following assumptions:Figure 7 shows the computational domain. The region represents the scaffold wall, the open interval is the growing tissue whereas the region corresponds to the interstitial fluid that brings nutrient to the growing construct. We denote by the boundary of the computational domain and by n the outward unit normal vector on . We have at and at .

Fig. 7

Schematic representation of the 1D scaffold-based bioreactor

all model variables depend on the sole spatial coordinate x and on the time variable t;

the solid displacement field has only one nonvanishing component, that is with ;

the strain tensor has only one nonvanishing component, that is ;

The 1D mechanobiological model consists of the following PDE subsystems:

Poroelastic IV-BVP for biomass for given and , , find the solid displacement , the fluid friction velocity and the pressure that satisfy the following system of partial differential equations in balance form: whereis the tissue permeability with defined in Eq. (10d) (right) whereas is the so-called aggregate modulus [69] with . To close the problem, we specify the following initial and boundary conditions:

Mass balance IV-BVP for nutrient concentration for given u, V and , , find the oxygen nutrient concentration that satisfies the following system of partial differential equations in balance form: where: and To close the problem, we specify the following initial and boundary conditions:

Mass conservation IV-BVP for cellular populations for given u, p and c, find the volume fractionsand that satisfy the following system of partial differential equations in balance form: where is the solid phase velocity, the production terms are: and the fluid fraction is computed asTo close the problem, we specify the following initial and boundary conditions: The boundary conditions (17i)- (17j) express the fact that cellular phases can flow out of the biomass only because of the presence of an advective field.

Indicator of the isotropy of the local stress state

In the 1D configuration, the total stress tensor can be decomposed into the sum of isotropic and anisotropic components given by: where is the unit vector . We use to measure the degree of anisotropicity of the stress state at any point x of the mixture and at any time t. Namely, we define the parameter r aswhere is the Frobenius norm. We can give a mechanical interpretation of (18c) by studying the Mohr circle at point (x, t). The principal components of are: from which it follows that the Mohr circle at (x, t) has center and radius equal to the maximum total shear stress at (x, t)Comparing the latter relation with (18c) we conclude that the indicator of the local stress state anisotropy can be written asWe also need to characterize an appropriate value for the threshold parameter representing the level of hydrodynamic shear stress that induces metabolic activity of the cell population and therefore separates the isotropic regime from the anisotropic regime. In [57, 58] it is shown that hydrodynamic shear below 10 mPa may promote GAG synthesis, so that, coherently with (18f), we assume

Numerical approximation of the 1D mechanobiological model

Solving in closed form the mechanobiological model proposed in this article is a very difficult task because of the strong nonlinear nature of the problem. Therefore, in this section we illustrate the approximation methods that are used to solve numerically the equation system in the 1D setting of Fig. 7.

Computational algorithm

Prior to discretization, we need to reduce the solution of the whole coupled system to the solution of a sequence of linearized equations of simpler form. For this purpose we set and subdivide the time interval into uniform subintervals of length , in such a way that the discrete time levels , , are obtained. Then, for each , we set and for all until convergence we perform the fixed point iteration schematically illustrated in Fig. 8. The algorithm consusts of two nested loops, a temporal outer loop and a spatial inner loop. For each step m of the inner loop we solve in sequence three linear PDE system blocks. For the convenience of the reader, for each substep of the inner loop ,we indicate in Fig. 8 the equations that are solved in the step referring to their numbering in Sect. 8. In order, the subproblems to be solved are: a poroelastic system for biomass displacement and fluid pressure, an advection–diffusion–reaction (ADR) system for oxygen concentration and an ADR system for cellular populations and ECM. For each system we indicate on the right of the corresponding block the values of the dependent variables that are given inputs whereas the updated values of the dependent variables that are returned as outputs of the block are indicated on the right of the downward arrow that exits out the block. We notice that as soon as a newly updated variable is available such variable is immediately plugged into the successive block as input variable. For this reason the algorithm of Fig. 8 can be regarded as a nonlinear block Gauss-Seidel method (see [52] Chapt. 7).

Fig. 8

Solution map

Two remarks are in order about the above described solution map. The first remark concerns the linear poroelastic system. The weak formulation of this problem leads to solving a saddle-point problem in block symmetric form to which the abstract analysis of [53], Chapt. 7 and [9] can be applied to prove existence and uniqueness of the solution pair . The second remark concerns with the two linear ADR problems to which the application of the maximum principle (see [59]) allows to prove nonnegativity of the solutions and , .

Finite element discretization

The computational procedure described in Sect. 9.1 leads to solving two kinds of BVPs: (1) a saddle-point problem; (2) two ADR equations. We numerically solve (1) and (2) using the Galerkin finite element approximation scheme on a family of partitions of the computational domain, h being the discretization parameter (see [53]). In the case of the saddle point problem (1) we employ piecewise linear finite elements on for both solid displacement and fluid pressure. Equal-order interpolation for u and p does not give rise to numerical instabilities as it would be the case if the Stokes equations for an incompressible fluid were to be solved (cf. [53], Chapt. 9), because in the present model the variable p is not a Lagrange multiplier (as in the Stokes system), rather, it is the solution of an elliptic Darcy problem (for a similar treatment see [12]). In the case of the ADR equation we employ for the approximation of the concentration and of the cellular volume fractions the primal-mixed finite element discretization scheme with exponential fitting stabilization proposed and investigated in [62]. This choice is taken because it ensures that the computed numerical solutions satisfy a strict positivity property even in the case of a strongly advective regime. Moreover, it can be checked that, if advective terms do not play a major role compared to oxygen molecular diffusion in the biomass, then the effect of the stabilization introduced by the primal-mixed method of [62] becomes negligible so that the accuracy of the scheme is not spoiled. This, instead, would not be the case if the classic upwind stabilization were adopted (see [10] for a discussion of this important issue).

Simulation tests

In this section we show the numerical results obtained by solving the 1D problem with the computational algorithm described in Sect. 9. We denote henceforth by and the normal stress and normal velocity externally applied at whereas is the external oxygen concentration. Two sets of simulation tests are performed. In the first set of simulations we set with the aim of investigating a static culture environment (see [57, 58]). In the second set of simulations we set and as in [14]. These values are characteristic of a culture in a perfusion bioreactor where an external hydrodynamic shear stress is applied [55, 58, 63].

The first investigated question is the effect of the input model parameter amount A of cell density at the beginning of the culture process () and at the pore wall (). For cells and ECM we set , , with , and for each set of simulations we use the following values of A: The above values of and agree with the biophysical evidence that at the beginning of the growth process, proliferating cells are present in larger amount than the other cellular populations.

, ;

, .

The second investigated question is the effect of the input model parameter cellular growth rate . In our computations we use two values of this parameter, and (cf. Table 1). These two values are selected by comparison with the maximum specific cell growth rate used in [63] in such a way that corresponds to “low growth regime” whereas corresponds to “high growth regime”.

The third investigated question is the effect of the input model parameter maximum value of the external oxygen concentration supplied to the growing structure by the surrounding environment. To determine the effect of oxygen availability on biomass growth we set and for all , and being the saturation and threshold oxygen concentration, respectively (see Table 1).

For a synthetic representation of the isotropy indicator r, we define the following (equivalent) parameter as if and if . In the remainder of the discussion, no plot is reported for the fluid volume fraction because this variable can be computed by post-processing using (8f). Simulations are run over the time interval , with days, and the one-dimensional plots show the time evolution of the solid and fluid mixture components at the spatial coordinate . The values of model parameters used in the numerical experiments are reported in Table 1.

Discussion of simulation results

In this section we address a critical discussion of the more significant outcomes of the simulations of the model illustrated in Sect. 8. The adoption of the 1D setting has three points of strength. The first point of strength is that the simplicity of the geometry permits a verification of reliability of model predictions based on biophysical intuition. The second point of strength is that, despite being simple, the 1D setting preserves the main features of the 3D biomass growth process, permitting a comparison with experimental measurements. The third point of strength is that it is relatively easy to single out the presence of critical parameters in the mathematical formulation and investigate their quantitative influence on the evolution of mixture components.

The illustrated numerical results indicate that the in vitro cell cultivation process is strongly sensitive to variations of (1) the initial seeding density of cells, (2) the value of the maximum growth rate and (3) the mechanical boundary conditions. In particular, the amount of seeded cells turns out to be determinant for cell responsiveness: initial cell density should be high enough to ensure optimal conditions for proliferation, but not so high that grow factors are rapidly depleted from the medium and the contact inhibition phenomenon prevents the formation of new colonies (see Fig. 13). To furtherly support this conclusion, Fig. 12 shows that at high seeding density no significant change in cell number is predicted by the model in the first stage of the culture because proliferating cells fluctuate around a mean value, that represents the average level of proliferation measured in the first 2 weeks of the experiments (see Fig. 3c of [11]). Furthermore the value of parameter influences the long-term behavior of the biomass. Model simulations indicate that if is smaller than the reference value , cell metabolic activity and ECM synthesis significantly decrease in the cultured construct and are completely exhausted at about 10 days of culture (Figs. 9, 11 top). On the contrary, model predictions show that if the maximum growth rate exceeds the reference value, cell and ECM volumetric fractions increase until convergence to a finite value that represents a stable steady state of the mathematical system (Figs. 10, 11 bottom). This finding represents a favorable result from the experimental point of view, because it predicts the formation, at the end of the cultivation and under specific conditions, of the bio-artificial texture to be used for replacing damaged tissues, that constitutes the real aim of TE. A similar objective can be reached by conveniently assigning the mechanical boundary conditions at the interface between the biomass construct and the interstitial fluid. Model results indicate that when the biomass is stimulated by both external fluid velocity and pressure, even if is tuned on a under-threshold value, the amount of cells and ECM in the construct remains considerable until the end of the simulation (see Fig. 14 top panel). Such a behavior is mainly due to the fact that the external force exerted by the fluid gives rise to an anisotropic stress state that instantaneously propagates throughout the domain, as shown in Fig. 16 bottom, maintaining the anisotropic mechanical configuration in both low and high growth regimes. This mechanical stimulus is the sole responsible of the strong mitotic functional activity occurring within the biomass because, as evidenced in Fig. 15, oxygen consumption is practically absent. This outcome reinforces the notion that mechanical stimulation in perfused cultures may promote chondrogenesis and ECM production [17, 57, 58]. Actually, in order to achieve this optimal result, nutrient concentration at the fluid-biomass interface should not fall under a critical level otherwise cell functionality could be rapidly reduced until cell apoptosis (see Fig. 17 top). In these conditions, cell survival is ensured only if the growth rate is large enough (Fig. 17 bottom).

Fig. 13

Spatial and temporal evolution of oxygen concentration in the static culture for . Initial condition IC2. Left: . Right: . In the case of initial condition IC1 we observe a similar behavior of except for a delayed decay

Fig. 12

Temporal evolution of cellular populations and ECM in the static culture for . Initial condition IC2. Top: . Bottom: . Solid line: ; dashed line: ; dotted line: ; dash-dot line:

Fig. 9

Temporal evolution of cellular populations and ECM in the static culture for . Initial condition IC1; . Solid line: ; dashed line: ; dotted line: ; dash-dot line:

Fig. 11

Temporal evolution of cellular populations and ECM in the static culture for . Initial condition IC1. Top: . Bottom: . Solid line: ; dashed line: ; dotted line: ; dash-dot line:

Fig. 10

Temporal evolution of cellular populations and ECM in the static culture for . Initial condition IC1. Top: . Bottom: , zoom of the first 8 days of culture. Solid line: ; dashed line: ; dotted line: ; dash-dot line:

Fig. 14

Temporal evolution of cellular populations and ECM in the perfused culture for . Initial condition IC1. Top: . Bottom: . Solid line: ; dashed line: ; dotted line: ; dash-dot line:

Fig. 16

Spatial and temporal evolution of parameter in the static culture (top) and in the perfused culture (bottom) for . Initial condition IC1. Top left: . Top Right: . Bottom:

Fig. 15

Spatial and temporal evolution of oxygen concentration in the perfused culture for . Initial condition IC1. Left: . Right:

Fig. 17

Temporal evolution of cellular populations and ECM in the perfused culture for . Initial condition IC1. Top: . Bottom: . Solid line: ; dashed line: ; dotted line: ; dash-dot line:

The behavior of the solid mixture components is in excellent agreement with the experimental trends obtained by cultivation of engineered tissues in bioreactors. In particular, the temporal evolution of construct cellularity and ECM content, especially for an under-threshold value of , agree with experimental results shown in several papers [17, 23, 48, 72].

The characterization of the (an)isotropicity of the biomass intrinsic stress state through the equivalent parameter demonstrated to be a successful strategy to model the mechanical regulation of culture progression and to link the mechanisms occurring at the micro-scale level to the macroscopic functioning of the growing tissue (see [43]). Model predictions indicate that the parameter is an effective indicator of the propagation of the isotropic and anisotropic waves within the construct and allows an easy and immediate identification of the adhesion mechanisms developing at the single cell-level, that, accordingly, drive the evolution of the volumetric fraction and respectively.

Conclusions and future perspectives

In the present article we propose a novel mathematical formulation based on the continuum assumption to describe the biomechanical sensitivity of articular chondrocytes. The natural application of our model is Tissue Engineering, a continuously growing discipline within the wider area of Regenerative Medicine, in which the control of cell response to multi-factorial stimuli is of utmost importance to obtain products suitable to clinical practice. However, it is worth noting that the proposed scheme may be used as well to describe more general settings in Cellular Biology, for example, the expansion of staminal cells.

The principal novelty of our contribution is the development of a model based on the use of Partial Differential Equations (PDEs) that incorporates the concept of “force isotropy” on the cell within the general and well established framework of poroelastic theory of mixtures and of cell population models. Specifically, the model translates into a simplified mathematical formalism, based on the use of suitably parametrized Heaviside functions, how the induced cytoskeletal tensional states trigger signalling transduction cascades regulating functional cell behavior, for example, the traslocation of specific transcription factors in the nucleus. According to the concept of force isotropy, it turns out that if cell adhesion-mediated traction forces have approximately the same strength over the cell surface, then the cell nucleus tends to maintain a roundish morphology, otherwise the cell nucleus tends to elongate. In the first case, the cell tensile condition is defined as “isotropic cytoskeletal tension” whereas in the second case the cell tensile condition is defined as “anisotropic cytoskeletal tension”.

Having defined the cytoskeletal stress characterization at the single cellular level, the next step of our approach is to build upon the concept of continuum-based approach to extend the above described description to the local stress tensor associated with the biological construct to mathematically represent the isotropic or anisotropic cell adhesion state. To this purpose, we generalize in a natural manner the previous definitions prescribing that if the anisotropic part of the local stress tensor is lower than a fixed threshold then the local stress state of the system is isotropic otherwise the local stress state of the system is anisotropic.

The final step of our model construction is to incorporate the above illustrated mechanobiological scheme within the setting of the theory of poroelasticity of a mixture composed by a solid and a multi-component fluid phases. The mixture represents the cellular construct in which several different cellular populations are well-mixed and oxygen delivery and consumption is taken into account to regulate in a dynamical manner the progressive fate of the evolving (macroscopic) tissue. The overall mathematical formulation consists of a system of conservation laws (mass and linear momentum) for the phases and components of the mixture that includes stress state and oxygen tension as main determinants of cellular culture evolution.

A thorough investigation of the PDE system is critically performed in a simplified 1D setting to allow an easy preliminary validation of the formulation. Extensive simulation tests outline a generally sound response of the computational model with respect to biophysical conjectures. In particular, numerical results indicate that the in vitro cell cultivation process is strongly sensitive to variations of (1) the initial seeding density of cells, (2) the value of the maximum growth rate and (3) the mechanical boundary conditions.

Below, we mention several future steps that we intend to take in the prosecution of this promising research activity.

A stability analysis of the homogeneous steady states of the dynamical system that describes the conservation of mass of the solid mixture components. Such a study will allow us to characterize the admissible range of values of model constitutive parameters that ensures the biophysical consistency of the proposed mathematical representation, in the same spirit as in [41] and [37].

The introduction of a visco-elastic component in the constitutive law for the total stress (as recently done in [7]). This extension of the model will allow us to perform a validation of the model and of the computational tool against available analytical solution and data (see [69]).

The inclusion of other mixture constituents, such as proteoglycan and collagen as done in [33]. This extension of the model will allow us to provide a more realistic biomechanical description of the growing tissue.

The extension of the computational algorithms to treat a fully three-dimensional representation of the scaffold pore to allow a deeper model validation against previous existing simulation results and experimental data (see, e.g., [19, 34, 56]). In particular, a 3D implementation of the model could provide a very interesting in silico scheme to simulate different levels of isotropy/anisotropy that otherwise should be reproduced in vitro at the price of complex engineering strategies, as described in [34]. Once a wide range of simulated mechanical configurations is available, it should be easier to relate cell response to the external mechanical stimuli and, consequently, to predict cell behavior in terms of cell adhesion, proliferation and differentiation. The 3D model could also be used to reproduce one pore of the microscaffold structures recently developed in [44] to mimic the native cellular environment. Cells confined into micropores are subject to similar environmental cues as in vivo, so that the behavior predicted by 3D computations could be directly compared with the in vivo cellular processes. Of course, passing to a 1D implementation to a fully 3D simulation tool requires to face and solve, at least, four computational challenges. The first challenge is the need of a flexible tool for the generation of an accurate geometrical description of the structure to be simulated. The preferable choice is to use tetrahedral elements and to this purpose a very good 3D mesh generator is the open-source program gmesh. The second challenge is the selection of a stable and accurate time-advancing discretization method. The obvious choice is to continue to emply the Backward Euler scheme. If a more accurate method is in order, the choice might fall on the second-order Trapezoidal or TR-BDF2 methods (for description and analysis, see [52], Chapter 11). The third challenge is the need of extending the selection of the stable and accurate finite element spaces used to approximate the various subproblems to be solved with the fixed-point map illustrated in Sect. 9.1. To this purpose, the Mini element and the Taylor-Hood pair are the best options for an accurate and stable treatment of the poroelastic equations, whereas the Edge Averaged Finite Element scheme investigated in [75] is a very effective method for accurately dealing with sharp layers in the nutrient concentration and/or cellular population profiles while ensuring the positivity of the computed solution. In the perspective of implementing the mechanobiological model within a 3D finite element framework a possible interesting programme might be to exploit the facilities of the software MP-FEMOS (Multi-Physics Finite Element Modeling Oriented Simulator) that has been developed by one of the authors [1, 39, 40, 61].