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.

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

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 neotissue. Cartilage tissue growth in engineered constructs had been already studied in a series of papers by Klisch and coauthors [30][31][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 surfacestress 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 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 factors which regulate cell behavior (we refer to [43] and references therein for a more biologically detailed discussion of these complex processes).
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).
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 massconserving 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: (a) 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]; (b) the poroelasticity theory is used to model the interaction of deformation and fluid flow in the fluid-saturated porous, elastic solid [6,22]; (c) 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][32][33]; (d) 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]; (e) 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 x and t denote the space and time variables, respectively. We use the convention that the dependence of all variables and model parameters on x and t is left understood except otherwise stated.
The geometrical configuration of the mixture is identified by the open bounded set X & R d (d ¼ 3 unless otherwise specified). The domain X does not evolve in time, rather, it is the amount of each species at a point x 2 X 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 Q T end :¼ X Â ð0; T end Þ the space-time cylinder in which the TE problem is studied, T end [ 0 being the final time of culture process.
Referring to Fig. 3, for all t [ 0 we associate with a generic point x 2 X a fixed representative elementary volume (REV) V x (see [73]) and denote by jV x j its ddimensional volume. Then, for each component i ¼ fl; s of the growing mixture, we define the volume fraction as 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 According to the biochemical hypothesis a), we consider three ACC populations: proliferating, ECM secreting state and quiescent state, denoted by the letters n, v and q, respectively. Then, we have Eventually, we denote by / ¼ / fl ; / n ; / v ; / q ; Â / ECM T 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.
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 q w 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 socalled 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 u s ¼ u s ðx; tÞ and v s ¼ The growth process of each mixture component (cellular growth and ECM secretion) is taken into account by introducing the following decomposition [33] e g ¼ e g g þ e e g ; g ¼ cells; ECM; ð2bÞ where e g g 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 e e g 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, v fl being the velocity of the interstitial fluid. For notational brevity, from now on, we simply write u and e instead of u s and e s , respectively. 4 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 Q T end : where J / 2 R 5Âd is the flux matrix and Q is the net mass production rate for which the following constraint holds, due to Assumption 3 X f¼s;fl 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 Q T end : T s ðu; p; /Þ ¼ X g¼cells;ECM / g T g ðu; p; /Þ ð4bÞ where r g is the effective stress tensor of the component g of the solid phase of the mixture, p ¼ pðx; tÞ is the pressure exerted by the fluid phase and I is the identity tensor. The isotropic stress ÀpI 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 T f , f ¼ s; fl, are the total stress tensors of the solid and fluid phases, while p f 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 p f we refer to [36] and [51]. We observe that, for all t 2 ð0; TÞ and at all x 2 X, it holds p s ðx; tÞ þ p fl ðx; tÞ ¼ 0:

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 v is the composite velocity of the mixture (cf. Eq. (2.4) of [51]). A simple manipulation allows us to write Eq. (5b) as In a similar manner, summing Eqs. (4) and using (4e), we get the total momentum equation Tðu; p; /Þ ¼ X where T 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) c ¼ cðx; tÞ 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 Q T end : the interstitial fluid velocity v fl being computed using (2c) as The mathematical description of the oxygen diffusion coefficient D c 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 coefficient where D c;fl and D c;s represent the nutrient diffusivity in the fluid and solid phase, respectively, while K eq 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 R g , g ¼ n; v; q, is the nutrient consumption rate of the cellular population / g and K 1=2 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 Q g , g ¼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 HðzÞ ¼ 0 for z\0 and HðzÞ ¼ 1 for z ¼ 0, z being a real number. Then, the stress-state dependent effect is represented by Hðr À rÞ, r and r being an indicator of the isotropy or anisotropy of the local stress state and a threshold value, respectively. If r\r, 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 Hðc À c thr Þ, c thr being a threshold concentration for cell activity. For notational brevity, we let H r :¼ Hðr À rÞ and H c :¼ Hðc À c thr Þ.
The following exchange/production rate terms are considered: -mitotic proliferation: we suppose cells in pool n proliferate at a rate given by Relation (7a) is a phenomenological law in which the term / n / fl keeps into account contact inhibition effects, while the term c K sat þ c k g is a nutrient-dependent modulation (Monod law [21]), K sat being the half-saturation constant and k g 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, V cell is the volume of a single cell, the constant of proportionality E [ 1 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 k GAG 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 / ECM reaches a maximum value / ECM;max 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 c thr , whereas apoptotic phase is related to age dependent cell death [60]. The time rates of change between state a (a ¼ n; v; q) and the inactive states (quiescence or apoptosis) are k qui and k apo , respectively; -exchange rate between pools n $ q: a first contribution in the direction n ! q is regulated by the mitotic characteristic inverse time constant 1=s m and take the form A second contribution is regulated by the probability rate b q!n that a cell in pool q enters into pool n, enhanced by the mechanical factor H r , giving the rate term -exchange rate between pools v $ q: the probability rate b v!q that a cell in pool v enters into pool q and the opposite for b q!n are mediated by the mechanical terms H r and 1 À H r , respectively, to signify that anisotropy favors proliferation while isotropy ECM secretion.
According to the exchange laws illustrated above, the production terms associated with cell populations are defined as: Q To conclude the mathematical description of mass exchange terms, we define the extracellular fluid production Q fl 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 / fl 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 postprocessing as 6 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 biomechanical 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 g h , g n 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 v, q and ECM components, while is accumulated along a specific polarization direction, identified by the unit vector d pol , 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 (v 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 (n 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 q and ECM are supposed to be isotropic. The definition of d pol in Eq. (9b) and the law for its time evolution is a delicate issue. In [8], d pol is characterized according to the dynamics of the evolution of the cell cytoskeleton, which reorganizes itself according to its mechanosensing mechanisms. A simplified version of the model proposed in [8], and also adopted in the present work, is represented by the choice where d e is the normalized eigenvector of the infinitesimal strain tensor e s , associated with the eigenvalue of largest module, which physically corresponds to the maximum principal dilatation of the biomass around such a point [68]. The coefficients g g ðx; t; /Þ, g ¼ n; v; q; ECM, give a measure of the amount of mass of the cellular population of type g deposited at time t at point x. 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 x 2 X and for g ¼ cells; ECM: The quantity c R;g ðx; t; /Þ represents the amount of mass of the cellular population g deposited at time t at point x 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 c R;g ðx; t; /Þ :¼ Q g ðx; t; /Þ; g ¼ cells; ECM, in such a way that the initial value problems that furnish the characterization of the growth coefficients become, for h ¼ v; q; ECM: g h ðx; t; /Þ ¼ 0 and o ot g n ðx; t; /Þ ¼ Q n ðx; t; /Þ t 2 ð0; T end ð9hÞ g n ðx; t; /Þ ¼ 0 ð9iÞ having used the identities TrðIÞ ¼ 3 and Trðd pol d pol Þ ¼ 1.

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 r g ðu; / g Þ ¼ 2l g e e g ðu; / g Þ þ k g Tre e g ðu; / g ÞI ¼ 2l g eðuÞ À e g g ð/ g Þ þ k g Tr eðuÞ À e g g ð/ g Þ I; where k g and l g are the Lamé parameters of each component of the solid phase, g ¼ n; q; v; ECM, and e g g ð/ g Þ 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 Kð/ fl Þ ¼ / 2 fl C F I is defined as in [65], C F being a friction coefficient. To provide a physically consistent characterization of C F we apply the classic Stokes theory for viscous drag to the biomass mixture and obtain R cell and l fl being cell radius and interstitial fluid dynamic viscosity, respectively, from which we get 7 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 X an open bounded set of R 3 representing the computational domain in which the biomass growth process physically takes place, and by C :¼ oX the boundary of X on which an outward unit normal vector n 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 U :¼ ðu; pÞ; c; / f gthe set of the dependent variables of the problem. Then, for each u 2 U we assume that C is divided into two disjoint portions, denoted by C u D and For a similar treatment of the splitting of the domain boundary we refer to [26] and [7]. 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 / ¼ /ðx; tÞ, find the solid displacement u, the fluid friction velocity w and pressure p such that the following equations are satisfied in Q T end : divTðu; p; /Þ ¼ 0 ð11bÞ Tðu; p; /Þ ¼ X g¼cells;ECM / g r g ðu; /Þ À pI ð11dÞ r g ðu; /Þ ¼ 2l g ðeðuÞ À e g g ð/ g ÞÞ þ k g TrðeðuÞ À e g g ð/ g ÞÞI supplied with the following initial and boundary conditions:  Mass balance IV-BVP for nutrient concentration: for given w ¼ wðx; tÞ, u ¼ uðx; tÞ and / ¼ /ðx; tÞ, find the nutrient concentration c and the nutrient flux density J c such that the following equations are satisfied in Q T end : supplied with the following initial and boundary conditions: Mass conservation IV-BVP for cellular populations: for given / ¼ /ðx; tÞ, c ¼ cðx; tÞ, T ¼ Tðx; tÞ and u ¼ uðx; tÞ, find the volume fractions / such that the following equations are satisfied in Q T end : supplied with the following initial and boundary conditions: where / 0 : X ! ðR þ Þ 5 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: -all model variables depend on the sole spatial coordinate x and on the time variable t; -the solid displacement field u has only one nonvanishing component, that is u ¼ ½u; 0; 0 T with u ¼ uðx; tÞ; -the strain tensor has only one nonvanishing component, that is e xx ðu; tÞ ¼ ouðx; tÞ=ox; Figure 7 shows the computational domain. The region x\0 represents the scaffold wall, the open interval X ¼ ð0; LÞ is the growing tissue whereas the region x [ L corresponds to the interstitial fluid that brings nutrient to the growing construct. We denote by C :¼ oX ¼ 0; L f g the boundary of the computational domain and by n the outward unit normal vector on oX. We have n ¼ À1 at x ¼ 0 and n ¼ þ1 at x ¼ L.
The 1D mechanobiological model consists of the following PDE subsystems: Poroelastic IV-BVP for biomass for given / g and g g , g ¼ cells; ECM; fl, find the solid displacement u : Q T end ! R, the fluid friction velocity w : Q T end ! R and the pressure p : Q T end ! R that satisfy the following system of partial differential equations in balance form: where is the tissue permeability with K ref defined in Eq. (10d) (right) whereas H A ¼ k þ 2l is the socalled 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 / g , g ¼ n; v; q; fl, find the oxygen nutrient concentration c : Q T end ! R þ that satisfies the following system of partial differential equations in balance form: where: and D c ¼ D c;fl 3k À 2/ fl ðk À 1Þ 3 þ / fl ðk À 1Þ ; To close the problem, we specify the following initial and boundary conditions: oc ox cðL; tÞ ¼ c ext ðtÞ: ð15iÞ Mass conservation IV-BVP for cellular populations: for given u, p and c, find the volume fractions and / fl : Q T end ! R þ that satisfy the following system of partial differential equations in balance form: where v s ¼ ou=ot is the solid phase velocity, the production terms are: Pð/; c; TÞ ¼ ð17cÞ Cðc; TÞ ¼ and the fluid fraction is computed as To 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 T can be decomposed into the sum of isotropic and anisotropic components T ¼ T iso þ T aniso given by: where d pol is the unit vector ½1; 0; 0 T . We use T aniso 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 as rðx; tÞ ¼ kT aniso ðx; tÞk F 2l ¼ / s ðx; tÞ ouðx; tÞ ox À g n ðx; tÞ/ n ðx; tÞ where k Á k F 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 T are: from which it follows that the Mohr circle at (x, t) has center C ¼ ðr I þ r II Þ=2 and radius equal to the maximum total shear stress at (x, t) s max ðx; tÞ ¼ r I ðx; tÞ À r II ðx; tÞ 2 ¼ l / s ðx; tÞ ouðx; tÞ ox À g n ðx; tÞ/ n ðx; tÞ : Comparing the latter relation with (18c) we conclude that the indicator of the local stress state anisotropy can be written as We also need to characterize an appropriate value for the threshold parameter r representing the level of hydrodynamic shear stress that induces metabolic activity of the cell population n 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 9 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 ½0; T end into N T ! 1 uniform subintervals of length Dt ¼ T end =N T , in such a way that the discrete time levels t k ¼ kDt, k ¼ 0; . . .; N T , are obtained. Then, for each k ¼ 0; . . .; N T À 1, we set U ð0Þ :¼ U n and for all m ! 0 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-diffusionreaction (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). 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 u ðmþ1Þ ; p ðmþ1Þ . 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 c ðmþ1Þ and / ðmþ1Þ g , g ¼ cells; ECM.

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 T h f g h [ 0 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 T h 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 T b and V b the normal stress and normal velocity externally applied at x ¼ L whereas c ext is the external oxygen concentration. Two sets of simulation tests are performed. In the first set of simulations we set T b ¼ V b ¼ 0 with the aim of investigating a static culture environment (see [57,58]). In the second set of simulations we set T b ¼ 100mPa and V b ¼ 50lms À1 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 (t ¼ 0) and at the pore wall (x ¼ 0). For cells and ECM we set / g ðx; 0Þ ¼ A g expðÀx=L d Þ, g ¼ n; v; q; ECM, with L d ¼ L=5, and for each set of simulations we use the following values of A: The above values of A n and A g 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 k g . In our computations we use two values of this parameter, k g1 and k g2 (cf. Table 1). These two values are selected by comparison with the maximum specific cell growth rate k g0 used in [63] in such a way that k g \k g0 corresponds to ''low growth regime'' whereas k g [ k g0 corresponds to ''high growth regime''.
The third investigated question is the effect of the input model parameter maximum value of the external oxygen concentration c ext supplied to the growing structure by the surrounding environment. To determine the effect of oxygen availability on biomass growth we set c ext ðtÞ ¼ c sat and c ext ðtÞ ¼ c thr for all t 2 ½0; T end , c sat and c thr 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 n ¼ nðrÞ as nðrÞ ¼ 1 if r r and In the remainder of the discussion, no plot is reported for the fluid volume fraction / fl because this variable can be computed by post-processing using (8f). Simulations are run over the time interval ½0; T end , with T end ¼ 30 days, and the one-dimensional plots show the time evolution of the solid and fluid mixture components at the spatial coordinate x ¼ L=2. 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.
1. 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 k g influences the long-term behavior of the biomass. Model simulations indicate that if k g is smaller than the reference value k g;0 , 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 k g is tuned on a under-threshold value, the amount of cells and ECM in the construct remains considerable until 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 k g is large enough (Fig. 17  bottom).
In the case of initial condition IC1 we observe a similar behavior of c ox except for a delayed decay 2. 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 k g ,   [17,23,48,72]. 3. The characterization of the (an)isotropicity of the biomass intrinsic stress state through the equivalent parameter n 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 n 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 / v and / n 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 adhesionmediated 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 continuumbased 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 multicomponent 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. 1. 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]. 2. 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]). 3. 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. 4. 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].