Computational Modelling of Tissue-Engineered Cartilage Constructs

Cartilage is a fundamental tissue to ensure proper motion between bones and damping of mechanical loads. This tissue often suffers damage and has limited healing capacity due to its avascularity. In order to replace surgery and replacement of joints by metal implants, tissue-engineered cartilage is seen as an attractive alternative. These tissues are obtained by seeding chondrocytes or mesenchymal stem cells in scaffolds and are given certain stimuli to improve the establishment of mechanical properties similar to the native cartilage. However, tissues with ideal mechanical properties were not obtained yet. Growth and remodelling (G&R) computational models of tissue-engineered cartilage are invaluable to interpret and predict the effects of experimental designs. The current model contribution in the field will be presented in this chapter, with a focus on the response to mechanical stimulation, and the development of fully coupled modelling approaches incorporating simultaneously solute transport and uptake, cell growth, production of extracellular matrix and remodelling of mechanical properties.


4
• Section 2 emphasizes the transport, uptake and production of relevant metabolites or growth factors that impact the biosynthetic activity of chondrocytes or mesenchymal stem cells, and how thee mechanisms are a ffected by external stimuli. A particular focus will be given to the main metabolites involved in chondrocyte metabolism: glucose, oxygen and lactate; • Section 3 is related to the different models proposed to modulate the proliferation, death and migration of chondrocytes and, in the case of MSCs, the proliferation of these and their differentiation into chondrocytes and other possible lineages. • Section 4 is concerned with to the models of synthesis of the main compo nents of the extracellular matrix (ECM), glycosaminoglycans (GAGs) and collagen taking into account the impact of different stimuli on the production rates, binding and degradation of the matrix, as well as the alignment of collagen fibres to establish the anisotropy of the cartilaginous tissue. • Section 5 describes the models of the mechanical behaviour of cartilage and the remodelling of the mechanical properties of tissue engineered cartilage based on the produced biomass and ECM. • Section 6 presents the models that couple all the aforementioned con cepts into simultaneous metabolic, biosynthetic and mechanical remodelling models.

Models for solute transport, uptake and release
In order to obtain tissue engineered cartilage with a sufficient amount of extracellular matrix, cells need to consume high amounts of nutrients to support their anabolic activity. However, there are serious limitations to nutrient transport across the tissues, which become hindered by the increase of matrix accumulation and decrease of the tissue porosity. This limitation is particularly seen in the cores of the tissues, where supply of nutrients is limited and accumulation of toxic byproducts leads to increased heterogeneities in the growing tissues. As a consequence, tissues with cores with lower cell viability and ECM content are formed, leading to inferior mechanical properties [52,65,5,2]. It has been postulated that, apart from the diffusive transport present in free swelling and unstrained constructs, advective nutrient transport may be helpful to reduce heterogeneities in nutrient supply, with a higher positive contribution for the transport of large solutes over small solutes [52,80,12].
The simplest modelling approach for solute transport is based on the diffusion-reaction equation, where solute diffuses through a porous tissue with a diffusivity that is a fraction of the diffusivity in the fluid phase, and with a reactive term correspondent to the consumption or release of solutes depending on the amount of solute and the cell density in the tissue. A typical way to demonstrate the decrease of diffusivity across a porous tissue is given by the Mackie-Mears relationship. The most common representation of the reactive term is based on the Michaelis -Menten kinetics, as shown below [65,80,81,82,22,70]. When dynamic loading is involved, an advective term is included to represent the fluid flow mediated transport [81,70,25,68,71]. coefficient of the nutrient in the tissue, R the flux of the metabolite, ρcell the cell density, Vmax the maximu m uptake rate and Km is the half maximu m rate concentration and nf the fluid volume fraction. The solutes consumed by the cells that are typically simulated in previous works are glucose and oxygen and several studies used this simplified assumption to consumption with good results [81,65,22,68,9]. However, particularly in cases where the culture medium has a high content in glucose, the deleterious effect of lactate production in cell proliferation cannot be ignored and the release of lactate to the culture media is also modelled, considering both anaerobic and aerobic degradation of glucose depending on the experimental lactate to glucose ratios [65,36]. Under dynamic loading conditions, deformation affects solute transport in several ways. In first place, the diffusion coefficient depends on the porosity of the scaffold used for cell support. As stated before, a common law used in previous modelling contributions is the Mackie-Mears diffusion law. Another effect of dynamic loading in the cellular metabolism is related to the variation of cell density. Assuming, as a simplification, that the number of cells in a given volume is constant, under loading, the cell density is affected due to the change of volume of the constructs in a compressible scaffold. This volume change is described by the determinant of the deformation gradient tensor (J). This value describes the ratio between the volume of the deformed configuration and the undeformed configuration. Therefore, the deformed cell density is obtained as such [39,49].

Models for cellular dynamics
The cell population in tissue engineered cartilage is highly dynamic and dependent on several metabolic and physical cues. The main mechanisms associated with cells populating the newly formed tissue are [67]:  Proliferation -A fundamental factor to obtain ECM in proper amounts. However, the higher the cell population, the more likely the nutrient depletion and inhomogeneity in cell distributions across the tissue.  Differentiation -When the tissue is seeded with mesenchymal stem cells (MSCs), an important factor to control besides their proliferation is the differentiation into chondrocytes. Since MSCs can also give rise to adipocytes and osteocytes, a precise control of the biochemical and biomechanical cues to favour differentiation into a given precise lineage is fundamental.  Migrationspatial redistribution of cells in the scaffold can both occur due to random walks without a preferential direction or can happen directionally towards chemoattractants. In high density scaffolds, cells can form colonies.  Death -Apart from the regular lifespan of chondrocytes, lack of nutrients or aggressive external physical cues will speed up the process of death.
The typical model for cell dynamics in a tissue engineered cartilage takes into account these factors as follows: is introduced due to the assumption that new chondrocytes have mobility due to random walks [18]. While the death rate per cell, Rdeath, is assumed constant, the proliferation rate, Rprol, is modulated both by metabolic and mechanical factors that decrease the actual proliferation from the maximu m proliferation rate, µmax. The simplest model for nutrient-limited cell proliferation (accounted for in Rprol) in tissue engineered cartilage is given by the Monod kinetics. In this model, similar to the Michaelis-Menten kinetics for nutrient dynamics, growth is limited by the availability of a nutrient, for which a half-rate concentration controls the steepness until maximu m growth.
Several models based on the Monod kinetics have good results in comparison with experimental data, either by using one solute only, such as glucose [30,14] or oxygen [20,62,50,28,45], or a combination of solutes, like models with glucose and collagen [17], or models inhibited by pH decrease simplified as accumulation of lactate [36].
Another model that is commonly used in the literature to describe the limitation of chondrocyte growth by substrate is the Contois kinetics. This representation differs from the Monod kinetics because the growth in the Contois kinetics is also inhibited by the cell density, implying saturation of growth due to spatial competition of cells for resources, as shown below: The Contois kinetics has also provided good agreement to the growth of chondrocytes in different scaffolds, using glucose [30,18,15,9], oxygen [19,47], a combination of glucose and lactate [79] and a combination of glucose and lactate accounting for pH negative effects [16]. In terms of the metabolic modulation of cell growth, other less used equations are reported, such as the Heaviside step function [46,43], the Moser or heterogeneous n-th order model [30] and logistic function [68]. In order to incorporate the impact of mechanical stimulation in cell growth, some models have built up from the aforementioned mechanical factors and introduced the impact of shear stress in cell growth. It was shown experimentally that articular chondrocytes show a dose and time dependent response to shear stress [44,60]. Two main modeling contributions have been proposed to incorporate this effect in mathematical models thus far. In first place, a simple linear model of a linearly increase in the growth rate with increasin g shear stress was proposed [62,36]. An extension of this equation was proposed as a polynomial dependence with a non -integer factor by [47]. A more recent contribution resides in a piecewise function with shear stress, accounting for a maximu m stimulatory range of shear stresses between 0.1 and 0.6 Pa as determined experimentally and assuming suppression of growth for stresses above 1 Pa [56]. All the models presented thus far consider growth on a homogeneous cell population with the same characteristics . Other models have focused on particular compartments of the cell population for modeling, with a cell in a given state having a different role in tissue homeostasis. A model reporting a proliferative, an extracellular matrix prod ucing and a quiescent cell fraction was proposed with interchangeability between these compartments [63], having been recently expanded to include a transitional state between proliferative and ECM producing states and the possibility of quiescence and apoptosis [45]. Another compartmental modeling approach is related to the influence of the phase of the cell cycle during mitosis on the maturity of the cell and the possibility of undergoing protein synthesis [63]. While these contributions are valuable and in closer agreement with the inherent biology of the chondrocytes, the previously reported general chondrocyte growth models have shown good agreement with experiments and are, in most cases, a reasonable modeling approach.

Models for ECM growth
The extracellular matrix of articular cartilage is a collagen fibre network, mainly composed by type II collagen, and of glycosaminoglycans (GAGs) that provide mechanical support to chondrocytes and resistance to the mechanical stimuli that cartilage is subjected to [74,1,33]. In tissue engineered cartilage, it is highly required to stimulate the new tissues to produce a network similar to the native cartilage, with a content of 5-10% of the total mass in GAGs and 10-20% of the total mass in collagen [1]. While there are reported studies able to produce GAGs in a concentration similar to the native cartilage, the collagen content is much lower than the native values, being this one of the most significant hurdles to surpass to obtain viable tissue for implantation [10,11].
The main mechanisms behind ECM dynamics are:  Synthesis -Synthesis models are focused on the total cell population as a whole, assuming that all cells have the ability to produce ECM, or on a ECM producing cell compartment, depending on the model for cell growth used. Synthesis rates may depend on the availab ility of a given substrate or on the mechanical stimuli that cells face.  Binding -The newly synthesized ECM composing molecules are released into the media and then linked to the ECM.  Degradation -Bound ECM molecules are due to degrade and diffuse into the culture media at a given rate due to several forms of damage.
The first models appearing in the literature on this matter started to consider the ECM as a whole and did not provide distinction between collagen and proteoglycans. These models considered a linear synthesis rate modulated by the difference between current ECM concentrations and the steady -state concentration of ECM, due to the experimental observation that synthesis rates decay with the accumulation of ECM [74]. This model formulation was later adapted to include the impact of cell density in the growth rates and a separation by type of ECM component, the differentiation between bound and unbound ECM, and the rates of degradation [64,33,55,9,24,23,3]. The typical formulation for these three ECM mechanisms is depicted in the equations below: In the equations above, the unbound (ECMub), bound (ECMb) and the degradation products (ECMd) are controlled by the respective diffusion coefficients and by the synthesis (kECM,s), binding (kECM,b) and degradation rates (kECM,d). It is assumed that, in tissue engineered cartilage, ECM growth will saturate at a given concentration, hence the dependence of the synthesis rate on the steady state concentration ECMb,ss. Apart from the commonly used linear dependence of the synthesis rate on the concentration of the specific ECM entity to model, other dependences were reported, such as a logistic dependence on the unbound GAG concentration [58,57,45] or the impact of levels of relevant solutes, like glucose [7,54] and oxygen [58,57,45,7], on the s ynthesis rates of GAG. For collagen, a dependence on the cell proliferation time derivative, instead of the typical linear dependence on cell density, was also reported [5].
The impact of mechanical stimuli on the synthesis rates of extracellular matrix has been less explored so far but some contributions are provided in the literature. Fluid velocity levels were considered in a previous level to directly affect the GAG synthesis rate, based on experimental observations that fluid velocity has a stimulatory effect on GA G synthesis, while for collagen an attempt in the same work was performed considering augmented synthesis when the maximu m principal strain is above a given threshold value [78]. The dependence on fluid velocity for GAG s ynthesis and on maximu m principal strain for collagen synthesis was implemented with a different formulation in another work as well [7,8]. Fluid velocity and shear stress were also considered as stimulatory for cartilage growth above a given threshold for both proteoglycans and collagen [27]. More recently, a modular function assuming that there is an optimal cell volume for synthesis of GAG by chondrocytes was implemented as a function of tissue deformation [31].
Another factor that is relevant for the establishment of extracellular matrix and of the anisotropic properties is the remodeling and reorientation of the collagen fibers. It is known that collagen fibers align in preferred directions between the maximu m principal strain directions [75,3,41], therefore the application of mechanical stimuli can be used to drive the desired orientation of tissue engineered cartilage. The first work that applied the remodeling theory of collagen fibers in cartilage assumed that collagen fibers rotated with an angular velocity controlled by the angle between the collagen fiber directions in the undeformed configuration and the preferred fibril directions, taking into account as well the magnitude of the three possible spatial principal strains to establish these directions [75]. Other approach to computational modeling of fibre reorientation was provided through an anisotrophy tensor describing the degree of structural anisotrophy and an ellipsoid representation for the fibre material parameters. Here the reorientation is described through the angle between the current anisotrophy tensor and the Cauchy -Green strain tensor and mediated through a time constant [55]. Furthermore, a probability density approach for the distribution of collagen fibers that can change over time given mechanical stimuli was also proposed [26]. These concepts are inherently coupled with the constitutive relationship used to describe the mechanical behavior of growing cartilage.

Models for description of the mechanical behaviour and remodelling of mechanical properties
Cartilaginous tissue obtained through support of a porous scaffold is a mate¬rial with a very high water content, like the native articular cartilage. For this reason, a simple monophasic constitutive material is not sufficient to expla in the viscoelastic behavior of cartilage, caused by fluid flow-dependent and independent mechanisms inherent to the properties of the solid material [37]. For this reason, mixture models based on the biphasic theory were proposed. These models describe the total stress in the tissue with a solid stress σs and a hydraulic pore pressure component p [37,59,4,7,8,9]. , of the porous scaffold, accounting as well for the porosity, n, of the material [64].
Due to the very high fluid content of these native cartilage, both the solid and fluid phases are generally described as incompressible, or nearly incompressible for simplicity (that is, with a Poisson ratio close to 0.5). However, in the tissue engineered cartilage, some polymers are described as compressive solids in equilibrium with an incompressible fluid [53,69,41,9]. These models can partly describe the fluid flow related viscoelasticity related to the low permeability of the material.
As an extension to the biphasic model for cartilage behavior, triphasic models were developed to take into account the swelling behavior due to gradients in osmotic pressure. This model accounts for a fluid with ionic particles, inducing or limiting chemical expansion of the negatively charged proteoglycan chains due to electrostatic repulsion. The osmotic pressure gradient (Δπ) and the chemical potential of the fluid (µf ) that drive the ionic phase stress contribution is caused by differences in ion concentrations of the cartilage and the surrounding fluid [37].
For the solid phase of these models, several constitutitve relationships were proposed to describe the mechanical behavior. The simples theory to describe the behavior of porous and viscoelastic materials is the poroelastic theory, where the solid phase is linear eliastic and the fluid is viscous. In this model, the stress -strain relationship of the solid phase is provided by the Hooke Law: In the equation above, Ha is the aggregate modulus, which is a measure of the stiffness of the material in equilibriu m when fluid flow through the material ceases. This quantity is related to the Young's modulus (E) and the Poisson's coefficient (v) through the following relationship: For mechanical modelling of scaffolds impregnated with chondrocytes, the poroelastic theory has been widely used [52,66,69,9]. Another reported theory for modelling of the solid phase is the porohyperelastic theory. Similarly to the poroelastic theory, the fluid is viscous but the solid phase has an hyperelastic constitutive relation ship, such as the stress-strain behavior is modelled by a non-linear relationship dependent on the strain energy (W ) and on the deformation gradient tensor (F ) and its determinant.
Models with a solid phase described by the neo-Hookean hyperelastic relationship have been widely quoted to modulate the mechanical behavior of articular cartilage [77,80,72,4]. Modeling of the hydrogel solid phase has used both the neo-Hookean law or the Odgen law [40,13].
The hyperelastic models referred to until now are isotrophic with the same mehcanical properties in all dimensions. However, articular cartilage is an anisotrophic material, with stress -strain behavior dependent on the orientation of the collagen fibres [37]. Therefore, it is more appropriate in long-term studies of cartilage growth to model the growing tissue with an anisotrophic model with augmented tensile response in the loading directions equal to the fibre directions. A proposed model for this is the Holzapfel-Gasser-Odgen model [59], which divided the ECM into a nonfibrilar component, explained by the neo-Hookean model, and a fibrilar component with strain energy dependent on material parameters and the degree of anisotrophy of the tissue. A similar relationship was reported, with the difference of taking into account a continuous exponential angular fibre distribution [26,55].
Another theory reported to describe the mechanical behavior of hydrogels is the poroviscoelastic model, with a non-viscous fluid and a viscoelastic solid phase (Kalyanam et al. 2009;Roberts et al. 2011). This model is not used, to our knowledge, to model hydrogels with growing cartilage.
A highly relevant parameter for the description of the biphasic behavior of cartilage is the hydraulic permeability . This parameter is related to changes in the porosity and void ratio of the material. Several exponential relationships between the permeability and porosity or void ratio were presented, with two of the most common ones being the Holmes & Mow law [71,32,35] and the Carman-Kozeny law [62,9] for isotropic permeability remodeling. However, with the growth of collagen fibers, the permeability also becomes anisotropic, with different values according to the orientation parallel or perpendicular to the fibers. Studies for modeling of articular cartilage explants have already included this dependency [59,26].
While the newly formed tissue is growing and ECM is deposited, the mechanical properties of the tissue are changing. The target average values of mechanical properties of tissue engineered cartilage are between 450 to 800 kPa of compressive Young's modulus and 10 -16 to 10 -15 m 4 .N -1 .s -1 in hydraulic permeability [51]. A tissue that combines these two ranges of mechanical parameters has not yet been established. In long -term tissue engineered cartilage modeling, establishing remodeling algorithms to simulate and account for the change of the mechanical properties is fundamental to determine with accuracy the intrinsic mechanical response of the tissues to external stimuli.
Few studies have reported relationships for the modeling of the solid matrix properties under linear elastic assumptions. The young's modulus remodeling was previously described by a linear model for the aggregate modulus with the concentration of GAGs and collagen derived from experimental data on bovine cartilage in d ifferent ages [73]. An extension of this relationship was reported as well with a 4th order polynomial dependence on the collagen concentration [9]. If the Poisson's ratio is assumed constant, the Young's modulus can be derived directly from such relationships. However, a possible remodeling relationship for the Poisson's ratio related to the porosity of the material was adapted from [76], since a compressible material, with the growth of ECM, tends to approach incompressibility.
In anisotropic models, the remodeling of the non-fibrilar part is controlled by the concentration of GAGs and the remodeling of fibrilar part controlled by the concentration of collagen. One reported relationship relates the rate of remodeling with the ratio of the current concentrations to the expected steady state concentrations [55]. Regarding the remodeling of permeability, the lower availability of experimental permeability measurements compared to the modulus measurements hinders the fitting to mechanistic models. However, as a proper simplification, previously reported models relate the decrease of permeability with the increase of the volumetric fraction of cells and ECM throughout the construct, leading to a decrease in porosity [5].

Coupled metabolic and mechanical remodelling models
Currently, most of the modeling contributions for tissue engineered cartilage are focused in up to three of the modeling dimensions presented. While all these contributions are very valuable, for a complete description of the behavior of the tissue and the time and spatial evolution of mechanical properties, all four dimensions need to be included. The creation of a validated model with explanatory and predictive power with the dimensions of solute t ransport and consumption, cell dynamics, extracellular matrix growth and remodeling of mechanical properties will allow to explain in a more quantitative way the histologically observed differences in the distribution of the modeled quantities across the tissues, as well as being invaluable to recommend changes to the processes of tissue culture in order to obtain better results. Finally, the complete coupled model can be used to predict the impact of envisioned changes to the culture protocol, such as the type and geometry of the scaffold material, dynamic loading, culture exchange, seeding densities, among others. The general full modeling scheme flow is represented in Figure 1.

Fig. 1
Modeling scheme workflow for fully coupled tissue engineered cartilage growth and remodeling. Scheme employed by [7,8,9]. This work was expanded to a long-term culture case by simulating an experiment with both a solid and a channeled 2% w/w agarose construct with chondrocytes during 56 days. The model was calibrated with solid construct data on GAG and collagen concentrations, as well as on the compressive Young's modulus and validated by reproducing well the experimental data for the channeled construct [9]. This modeling effort allowed to gain quantitative insights on the spatial heterogeneity of the constructs, showing that the degree of spatial heterogeneity of the Young's modulus the constructs with a central diffusion channel is 23% of the control value, while for permeability the heterogeneity is 27% of the control one, showing a significant improvement for the channeled condition ( Figure 3). The degree of spatial heterogeneity and the insufficient permeability remodeling in the simulated solid constructs affects significantly the mechanical response to compressive strain, with nominal stresses for the simulated heterogeneous TE cartilage 57% lower than for native articular cartilage and pore pressures 53% lower than the native case [9]. Therefore, permeability is the main parameter to be improved to get a more similar mechanical response, c alling for new scaffold material designs and stimulation protocols. The developed complete model was also applied to simulate dynamic loading conditions. A model parametrized with literature parameters was used to simulate the distribution of cell density and ECM in cubic constructs subjected to either compression, shear or bending at 5% of height, 1Hz for 6h continuously [7]. While the simulation time is very short for relevant differences in the mechanical properties to be seen, bending was, under these conditions, the more favourable regime for cell proliferation and the spatial distributions are relevant with the establishment of cart ilage with different structural organizations due to different maximu m principal strain directions ( Figure 4). Current work is related to the model validation for the estimation of cell proliferation, ECM growth and mechanical properties remodelling under several different regimes of cyclic unconfined compressive loading. All the reported works until now have been focused on articular cartilage. However, the complete modeling scheme was also successfully applied to estimate the remodelling of temporomandibu lar joint disc, an area where tissue engineering is still in a very early phase. The application of static hydrostatic pressure for 72h on PEDGA -condylar constructs promoted a very slight improvement of the mechanical properties. This preliminary study pro vided a future basis to estimate the impact of long term static or dynamic hydrostatic pressure on the growth of new temporomandibular joint discs [8] The presented model applications for simulation of growth and remodeling under dynamic loading represent shortterm loading with limited differences in the mechanical properties between conditions. Future work involves the simulation of long-term intermitted compressive loading. On long-term regimes, the collagen content is relevant in terms of fibre organization, therefore the expansion of the constitutive relatio nships to include the anisotropic behaviour of tissue engineered cartilage and its time dynamics is a future goal.

Concluding remarks
This chapter provided an overview of the current state of the computational models used for simulation of growth of tissue engineered cartilage, namely the models based in mixture theory. The several individual contributions for the underlying phenomena, s uch as solute transport and uptake, cell growth, production of extracellular matrix and remodelling of mechanical properties were presented, with a focus on, when applicable, the models that incorporate the impact of mechanical stimulation on these phenomena. Several constitutive relationships to model the mechanical behaviour of tissue engineered cartilage have been proposed and integrated into these models. A fully coupled modeling approach was developed to accommodate all these phenomena in a simultaneou s fashion for a more realistic representation of the biomechanical phenomena and used to estimate the growth and remodelling of mechanical properties under free swelling and mechanical loading of tissues. Future challenges on this area include refinement o f the developed equations through validation with experimental studies with proper measurements of all the underlying variables when possible, study the degree of robustness and/or specificity with different cell-matrix systems, simulation of long term mechanical loading cultures and the incorporation of anisotropic models with reorientation of collagen fibers in the fully coupled formulation. 19 Articular cartilage is a fundamental tissue that resides in the surface of bones, providing a smooth and lubricated surface for relative bone motion in the joints and for transmission of loads with low friction [29,1]. Articular cartilage is generally between 2 and 4 mm thick. Unlike other tissues , it does not have surrounding blood vessels or nerves [29]. The cartilage is populated by chondrocytes, specialized cells for the production of extracellular matrix (ECM). This matrix is mostly composed of collagen fibers, proteoglycans, water and other less present components, such as noncollagenous proteins and glycoproteins. The ECM components are fundamental for water reten tion in the tissue, promoting a softer load transfer and motion [1,29]. Water is the most abundant component of articular cartilag e, generally accounting for 65 to 80% of the total tissue weight. The collagen of the articular cartilage is mostly type II and is related to the tensile resistance of cartilage. Collagen accounts for 10-20% of the total cartilage mass. The proteoglycans are composed of a protein core with glycosaminoglycans (GAGs) attached, being 5-10% of the cartilage mass. The GAGs have a global negative charge that helps to control the hydration of cartilage and to provide resistance to compression and shear. When under a mechanical load, the proteoglycans then promote a redistribution of water in cartilage, leading to an increase in osmotic pressure with water flow. The osmotic pressure becomes larger than the applied load, which is fundamental to protect bones from loading [1,29]. A fundamental characteristic of articular cartilage is the depth dependent organization, with three zones with distinct functions and collagen architectures: the superficial zone, responsible for protection against shear stress and with collagen fibers parallel to the surface of the tissue; the middle zone, with oblique collagen fibers and providing first resistance to compression; and the deep zone, with the highest resistance to compressive forces with collagen fibers oriented perpendicularly to the cartilage surface. The cartilage is anchored to the subchondral bone by the calcified layer [1].
The most common pathology associated with articular cartilage is os teoarthritis (OA), a degenerative disease that causes loss of the smooth surface of cartilage with pain, inflammation and loss of motion amplitude. The highest risk factor for OA is increasing age, while other factors such as obesity, genetics and gender are also associated [48]. The worldwide prevalence of OA was estimated to be 3.8% in 2010 [21] and the direct and indirect costs of the disease are very high. In the United States only, the annual medical care expenditures with OA are of about $185 billion [34]. While traditionally seen as a disease of the cartilage only, more recently OA has been identified as a multi-organ pathology, causing subsequent damage in bone marrow and bone, tendons, ligaments, muscles and neural tissues [29].
Since cartilage is an avascular tissue, the intrinsic regeneration capacity of articular cartilage is ve ry limited, leading to increasing severity of damage. The current therapeutic solutions are the total joint replacement by a metal implant , which is more common in older patients with very advanced damage. Other solutions for younger patients are the microfracture and autologous chondrocyte implantation to promote formation of new cartilaginous tissue in the injury site. These solutions have moderate short term success rates, while long term results are not satisfactory. The failure of these therapies is related with the formation of tissue with inferior mechanical properties to the native tissues, with possible fibrocartilage formation [42].
Tissue-engineered (TE) cartilage has been proposed as a prospective new treatment for osteoarthritis by the in vitro production of cartilaginous tis sue with more similar structure, composition and properties to the native articular cartilage. TE cartilage is obtained by seeding chondrocytes, or mes enchymal stem cells (MSCs) with chondrogenic cues, on a porous and biocompatible scaffold that is able to provide a favourable environment to maintain the differentiated phenotype of chondrocytes and to enable the production of extracellular matrix (ECM). Although promising, the translation of this approach to products has been hindered by factors such as insufficient mechanical properties, mainly due to the inability of the engineered tissues to have a type II collagen content similar to the native cartilage, difficulty in creating an anisotropic tissue structure with three layers with collagen fibres oriented as found in the native tissue, and heterogeneous mechanical properties with stiffer peripheries and softer cores [53,10,41,5].
In order to better predict the experimental conditions to subject the growing tissue to, either by mechanical, electrical, and chemical stimuli, computational models of tissue engineered cartilage are invaluable. Mathematical modeling in the context of TE cartilage has provided good insights on the nutrient distribution in the growing tissues [65,15,62,56], cell proliferation and death [15,62,56], synthesis of the main components of ECM, such as proteoglycans and collagen [74,33,55] and remodelling of biphasic mechan ical properties [73,15,62,55]. Most of these models attempt to solve one or two variables responsible for the full remodelling of TE cartilage. Recently, a new approach that couples all these factors in order to simulate spatiotempo ral patterns of metabolic activity, biomass growth and remodelling properties simultaneously was developed with results for both unloaded and mechanical stimulated constructs [7,8,9].
This chapter aims to review the body of work in the computational mod elling of tissue engineered cartilage with a focus on metabolic, biomass growth and mechanical remodelling. It is organized into several sections that emphasize different relevant aspects of the biomechanical behaviour of the growing cartilaginous tissues: • Section 2 emphasizes the transport, uptake and production of relevant metabolites or growth factors that impact the biosynthetic activity of chondrocytes or mesenchymal stem cells, and how thee mechanisms are a ffected by external stimuli. A particular focus will be given to the main metabolites involved in chondrocyte metabolism: glucose, oxygen and lactate; • Section 3 is related to the different models proposed to modulate the proliferation, death and migration of chondrocytes and, in the case of MSCs, the proliferation of these and their differentiation into chondrocytes and other possible lineages. • Section 4 is concerned with to the models of synthesis of the main compo nents of the extracellular matrix (ECM), glycosaminoglycans (GAGs) and collagen taking into account the impact of different stimuli on the production rates, binding and degradation of the matrix, as well as the alignment of collagen fibres to establish the anisotropy of the cartilaginous tissue. • Section 5 describes the models of the mechanical behaviour of cartilage and the remodelling of the mechanical properties of tissue engineered cartilage based on the produced biomass and ECM. • Section 6 presents the models that couple all the aforementioned con cepts into simultaneous metabolic, biosynthetic and mechanical remodelling models.

Models for solute transport, uptake and release
In order to obtain tissue engineered cartilage with a sufficient amount of extracellular matrix, cells need to consume high amounts of nutrients to support their anabolic activity. However, there are serious limitations to nutrient transport across the tissues, which become hindered by the increase of matrix accumulation and decrease of the tissue porosity. This limitation is particularly seen in the cores of the tissues, where supply of nutrients is limited and accumulation of toxic byproducts leads to increased heterogeneities in the growing tissues. As a consequence, tissues with cores with lower cell viability and ECM content are formed, leading to inferior mechanical properties [52,65,5,2]. It has been postulated that, apart from the diffusive transport present in free swelling and unstrained constructs, advective nutrient transport may be helpful to reduce heterogeneities in nutrient supply, with a higher positive contribution for the transport of large solutes over small solutes [52,80,12].
The simplest modelling approach for solute transport is based on the diffusion-reaction equation, where solute diffuses through a porous tissue with a diffusivity that is a fraction of the diffusivity in the fluid phase, and with a reactive term correspondent to the consumption or release of solutes depending on the amount of solute and the cell density in the tissue. A typical way to demonstrate the decrease of diffusivity across a porous tissue is given by the Mackie-Mears relationship. The most common representation of the reactive term is based on the Michaelis -Menten kinetics, as shown below [65,80,81,82,22,70].
When dynamic loading is involved, an advective term is included to represent the fluid flow mediated transport [81,70,25,68,71].
In the equations above, c represents the concentration of the nutrient, Dtissue represents the diffusion coefficient of the nutrient in the tissue, R the flux of the metabolite, ρcell the cell density, Vmax the maximu m uptake rate and Km is the half maximu m rate concentration and nf the fluid volume fraction. The solutes consumed by the cells that are typically simulated in previous works are glucose and oxygen and several studies used this simplified assumption to consumption with good results [81,65,22,68,9]. However, particularly in cases where the culture medium has a high content in glucose, the deleterious effect of lactate production in cell proliferation cannot be ignored and the release of lactate to the culture media is also modelled, considering both anaerobic and aerobic degradation of glucose depending on the experimental lactate to glucose ratios [65,36]. Under dynamic loading conditions, deformation affects solute transport in several ways. In first place, the diffusion coefficient depends on the porosity of the scaffold used for cell support. As stated before, a common law used in previous modelling contributions is the Mackie-Mears diffusion law. Another effect of dynamic loading in the cellular metabolism is related to the variation of cell density. Assuming, as a simplification, that the number of cells in a given volume is constant, under loading, the cell density is affected due to the change of volume of the constructs in a compressible scaffold. This volume change is described by the determinant of the deformation gradient tensor (J). This value describes the ratio between the volume of the deformed configuration and the undeformed configuration. Therefore, the deformed cell density is obtained as such [39,49].

Models for cellular dynamics
The cell population in tissue engineered cartilage is highly dynamic and dependent on several metabolic and physical cues. The main mechanisms associated with cells populating the newly formed tissue are [67]:  Proliferation -A fundamental factor to obtain ECM in proper amounts. However, the higher the cell population, the more likely the nutrient depletion and inhomogeneity in cell distributions across the tissue.
 Differentiation -When the tissue is seeded with mesenchymal stem cells (MSCs), an important factor to control besides their proliferation is the differentiation into chondrocytes. Since MSCs can also give rise to adipocytes and osteocytes, a precise control of the biochemical and biomechanical cues to favour differentiation into a given precise lineage is fundamental.
 Migrationspatial redistribution of cells in the scaffold can both occur due to random walks without a preferential direction or can happen directionally towards chemoattractants. In high density scaffolds, cells can form colonies.  Death -Apart from the regular lifespan of chondrocytes, lack of nutrients or aggressive external physical cues will speed up the process of death.
The typical model for cell dynamics in a tissue engineered cartilage takes into account these factors as follows: is introduced due to the assumption that new chondrocytes have mobility due to random walks [18]. While the death rate per cell, Rdeath, is assumed constant, the proliferation rate, Rprol, is modulated both by metabolic and mechanical factors that decrease the actual proliferation from the maximu m proliferation rate, µmax. The simplest model for nutrient-limited cell proliferation (accounted for in Rprol) in tissue engineered cartilage is given by the Monod kinetics. In this model, similar to the Michaelis-Menten kinetics for nutrient dynamics, growth is limited by the availability of a nutrient, for which a half-rate concentration controls the steepness until maximu m growth.
Several models based on the Monod kinetics have good results in comparison with experimental data, either by using one solute only, such as glucose [30,14] or oxygen [20,62,50,28,45], or a combination of solutes, like models with glucose and collagen [17], or models inhibited by pH decrease simplified as accumulation of lactate [36].
Another model that is commonly used in the literature to describe the limitation of chondrocyte growth by substrate is the Contois kinetics. This representation differs from the Monod kinetics because the growth in the Contois kinetics is also inhibited by the cell density, implying saturation of growth due to spatial competition of cells for resources, as shown below: The Contois kinetics has also provided good agreement to the growth of chondrocytes in different scaffolds, using glucose [30,18,15,9], oxygen [19,47], a combination of glucose and lactate [79] and a combination of glucose and lactate accounting for pH negative effects [16]. In terms of the metabolic modulation of cell growth, other less used equations are reported, such as the Heaviside step function [46,43], the Moser or heterogeneous n-th order model [30] and logistic function [68]. In order to incorporate the impact of mechanical stimulation in cell growth, some models have built up from the aforementioned mechanical factors and introduced the impact of shear stress in cell growth. It was shown experimentally that articular chondrocytes show a dose and time dependent response to shear stress [44,60]. Two main modeling contributions have been proposed to incorporate this effect in mathematical models thus far. In first place, a simple linear model of a linearly increase in the growth rate with increasin g shear stress was proposed [62,36]. An extension of this equation was proposed as a polynomial dependence with a non -integer factor by [47]. A more recent contribution resides in a piecewise function with shear stress, accounting for a maximu m stimulatory range of shear stresses between 0.1 and 0.6 Pa as determined experimentally and assuming suppression of growth for stresses above 1 Pa [56]. All the models presented thus far consider growth on a homogeneous cell population with the same characteristics . Other models have focused on particular compartments of the cell population for modeling, with a cell in a given state having a different role in tissue homeostasis. A model reporting a proliferative, an extracellular matrix prod ucing and a quiescent cell fraction was proposed with interchangeability between these compartments [63], having been recently expanded to include a transitional state between proliferative and ECM producing states and the possibility of quiescence and apoptosis [45]. Another compartmental modeling approach is related to the influence of the phase of the cell cycle during mitosis on the maturity of the cell and the possibility of undergoing protein synthesis [63]. While these contributions are valuable and in closer agreement with the inherent biology of the chondrocytes, the previously reported general chondrocyte growth models have shown good agreement with experiments and are, in most cases, a reasonable modeling approach.

Models for ECM growth
The extracellular matrix of articular cartilage is a collagen fibre network, mainly composed by type II collagen, and of glycosaminoglycans (GAGs) that provide mechanical support to chondrocytes and resistance to the mechanical stimuli that cartilage is subjected to [74,1,33]. In tissue engineered cartilage, it is highly required to stimulate the new tissues to produce a network similar to the native cartilage, with a content of 5-10% of the total mass in GAGs and 10-20% of the total mass in collagen [1]. While there are reported studies able to produce GAGs in a concentration similar to the native cartilage, the collagen content is much lower than the native values, being this one of the most significant hurdles to surpass to obtain viable tissue for implantation [10,11].
The main mechanisms behind ECM dynamics are:  Synthesis -Synthesis models are focused on the total cell population as a whole, assuming that all cells have the ability to produce ECM, or on a ECM producing cell compartment, depending on the model for cell growth used. Synthesis rates may depend on the availab ility of a given substrate or on the mechanical stimuli that cells face.  Binding -The newly synthesized ECM composing molecules are released into the media and then linked to the ECM.  Degradation -Bound ECM molecules are due to degrade and diffuse into the culture media at a given rate due to several forms of damage.
The first models appearing in the literature on this matter started to consider the ECM as a whole and did not provide distinction between collagen and proteoglycans. These models considered a linear synthesis rate modulated by the difference between current ECM concentrations and the steady -state concentration of ECM, due to the experimental observation that synthesis rates decay with the accumulation of ECM [74]. This model formulation was later adapted to include the impact of cell density in the growth rates and a separation by type of ECM component, the differentiation between bound and unbound ECM, and the rates of degradation [64,33,55,9,24,23,3]. The typical formulation for these three ECM mechanisms is depicted in the equations below: dependence on the cell proliferation time derivative, instead of the typical linear dependence on cell density, was also reported [5].
The impact of mechanical stimuli on the synthesis rates of extracellular matrix has been less explored so far but some contributions are provided in the literature. Fluid velocity levels were considered in a previous level to directly affect the GAG synthesis rate, based on experimental observations that fluid velocity has a stimulatory effect on GA G synthesis, while for collagen an attempt in the same work was performed considering augmented synthesis when the maximu m principal strain is above a given threshold value [78]. The dependence on fluid velocity for GAG s ynthesis and on maximu m principal strain for collagen synthesis was implemented with a different formulation in another work as well [7,8]. Fluid velocity and shear stress were also considered as stimulatory for cartilage growth above a given threshold for both proteoglycans and collagen [27]. More recently, a modular function assuming that there is an optimal cell volume for synthesis of GAG by chondrocytes was implemented as a function of tissue deformation [31].
Another factor that is relevant for the establishment of extracellular matrix and of the anisotropic properties is the remodeling and reorientation of the collagen fibers. It is known that collagen fibers align in preferred directions between the maximu m principal strain directions [75,3,41], therefore the application of mechanical stimuli can be used to drive the desired orientation of tissue engineered cartilage. The first work that applied the remodeling theory of collagen fibers in cartilage assumed that collagen fibers rotated with an angular velocity controlled by the angle between the collagen fiber directions in the undeformed configuration and the preferred fibril directions, taking into account as well the magnitude of the three possible spatial principal strains to establish these directions [75]. Other approach to computational modeling of fibre reorientation was provided through an anisotrophy tensor describing the degree of structural anisotrophy and an ellipsoid representation for the fibre material parameters. Here the reorientation is described through the angle between the current anisotrophy tensor and the Cauchy -Green strain tensor and mediated through a time constant [55]. Furthermore, a probability density approach for the distribution of collagen fibers that can change over time given mechanical stimuli was also proposed [26]. These concepts are inherently coupled with the constitutive relationship used to describe the mechanical behavior of growing cartilage.

Models for description of the mechanical behaviour and remodelling of mechanical properties
Cartilaginous tissue obtained through support of a porous scaffold is a mate¬rial with a very high water content, like the native articular cartilage. For this reason, a simple monophasic constitutive material is not sufficient to expla in the viscoelastic behavior of cartilage, caused by fluid flow-dependent and independent mechanisms inherent to the properties of the solid material [37]. For this reason, mixture models based on the biphasic theory were proposed. These models describe the total stress in the tissue with a solid stress σs and a hydraulic pore pressure component p [37,59,4,7,8,9]. , of the porous scaffold, accounting as well for the porosity, n, of the material [64].
Due to the very high fluid content of these native cartilage, both the solid and fluid phases are generally described as incompressible, or nearly incompressible for simplicity (that is, with a Poisson ratio close to 0.5). However, in the tissue engineered cartilage, some polymers are described as compressive solids in equilibrium with an incompressible fluid [53,69,41,9]. These models can partly describe the fluid flow related viscoelasticity related to the low permeability of the material.
As an extension to the biphasic model for cartilage behavior, triphasic models were developed to take into account the swelling behavior due to gradients in osmotic pressure. This model accounts for a fluid with ionic particles, inducing or limiting chemical expansion of the negatively charged proteoglycan chains due to electrostatic repulsion.
The osmotic pressure gradient (Δπ) and the chemical potential of the fluid (µf ) that drive the ionic phase stress contribution is caused by differences in ion concentrations of the cartilage and the surrounding fluid [37]. σ = σ s − (Δπ +µf )I (14) For the solid phase of these models, several constitutitve relationships were proposed to describe the mechanical behavior. The simples theory to describe the behavior of porous and viscoelastic materials is the poroelastic theory, where the solid phase is linear eliastic and the fluid is viscous. In this model, the stress -strain relationship of the solid phase is provided by the Hooke Law: In the equation above, Ha is the aggregate modulus, which is a measure of the stiffness of the material in equilibriu m when fluid flow through the material ceases. This quantity is related to the Young's modulus (E) and the Poisson's coefficient (v) through the following relationship: For mechanical modelling of scaffolds impregnated with chondrocytes, the poroelastic theory has been widely used [52,66,69,9]. Another reported theory for modelling of the solid phase is the porohyperelastic theory. Similarly to the poroelastic theory, the fluid is viscous but the solid phase has an hyperelastic constitutive relation ship, such as the stress-strain behavior is modelled by a non-linear relationship dependent on the strain energy (W ) and on the deformation gradient tensor (F ) and its determinant.
Models with a solid phase described by the neo-Hookean hyperelastic relationship have been widely quoted to modulate the mechanical behavior of articular cartilage [77,80,72,4]. Modeling of the hydrogel solid phase has used both the neo-Hookean law or the Odgen law [40,13].
The hyperelastic models referred to until now are isotrophic with the same mehcanical properties in all dimensions. However, articular cartilage is an anisotrophic material, with stress -strain behavior dependent on the orientation of the collagen fibres [37]. Therefore, it is more appropriate in long-term studies of cartilage growth to model the growing tissue with an anisotrophic model with augmented tensile response in the loading directions equal to the fibre directions. A proposed model for this is the Holzapfel-Gasser-Odgen model [59], which divided the ECM into a nonfibrilar component, explained by the neo-Hookean model, and a fibrilar component with strain energy dependent on material parameters and the degree of anisotrophy of the tissue. A similar relationship was reported, with the difference of taking into account a continuous exponential angular fibre distribution [26,55].
Another theory reported to describe the mechanical behavior of hydrogels is the poroviscoelastic model, with a non-viscous fluid and a viscoelastic solid phase (Kalyanam et al. 2009;Roberts et al. 2011). This model is not used, to our knowledge, to model hydrogels with growing cartilage.
A highly relevant parameter for the description of the biphasic behavior of cartilage is the hydraulic permeability . This parameter is related to changes in the porosity and void ratio of the material. Several exponential relationships between the permeability and porosity or void ratio were presented, with two of the most common ones being the Holmes & Mow law [71,32,35] and the Carman-Kozeny law [62,9] for isotropic permeability remodeling. However, with the growth of collagen fibers, the permeability also becomes anisotropic, with different values according to the orientation parallel or perpendicular to the fibers. Studies for modeling of articular cartilage explants have already included this dependency [59,26].
While the newly formed tissue is growing and ECM is deposited, the mechanical properties of the tissue are changing. The target average values of mechanical properties of tissue engineered cartilage are between 450 to 800 kPa of compressive Young's modulus and 10 -16 to 10 -15 m 4 .N -1 .s -1 in hydraulic permeability [51]. A tissue that combines these two ranges of mechanical parameters has not yet been established. In long -term tissue engineered cartilage modeling, establishing remodeling algorithms to simulate and account for the change of the mechanical properties is fundamental to determine with accuracy the intrinsic mechanical response of the tissues to external stimuli.
Few studies have reported relationships for the modeling of the solid matrix properties under linear elastic assumptions. The young's modulus remodeling was previously described by a linear model for the aggregate modulus with the concentration of GAGs and collagen derived from experimental data on bovine cartilage in d ifferent ages [73]. An extension of this relationship was reported as well with a 4th order polynomial dependence on the collagen concentration [9]. If the Poisson's ratio is assumed constant, the Young's modulus can be derived directly from such relationships. However, a possible remodeling relationship for the Poisson's ratio related to the porosity of the material was adapted from [76], since a compressible material, with the growth of ECM, tends to approach incompressibility.
In anisotropic models, the remodeling of the non-fibrilar part is controlled by the concentration of GAGs and the remodeling of fibrilar part controlled by the concentration of collagen. One reported relationship relates the rate of remodeling with the ratio of the current concentrations to the expected steady state concentrations [55]. Regarding the remodeling of permeability, the lower availability of experimental permeability measurements compared to the modulus measurements hinders the fitting to mechanistic models. However, as a proper simplification, previously reported models relate the decrease of permeability with the increase of the volumetric fraction of cells and ECM throughout the construct, leading to a decrease in porosity [5].

Coupled metabolic and mechanical remodelling models
Currently, most of the modeling contributions for tissue engineered cartilage are focused in up to three of the modeling dimensions presented. While all these contributions are very valuable, for a complete description of the behavior of the tissue and the time and spatial evolution of mechanical properties, all four dimensions need to be included. The creation of a validated model with explanatory and predictive power with the dimensions of solute t ransport and consumption, cell dynamics, extracellular matrix growth and remodeling of mechanical properties will allow to explain in a more quantitative way the histologically observed differences in the distribution of the modeled quantities across the tissues, as well as being invaluable to recommend changes to the processes of tissue culture in order to obtain better results. Finally, the complete coupled model can be used to predict the impact of envisioned changes to the culture protocol, such as the type and geometry of the scaffold material, dynamic loading, culture exchange, seeding densities, among others. The general full modeling scheme flow is represented in Figure 1.
For the simulation of free swelling constructs with different scaffold geometries, the model was applied to simulate short-term effects in the Young's modulus and hydraulic permeability of constructs with cylindrical and cubic geometries, either solid or with a central channel [6]. Despite the short culture period of 72h, it was possible to determine that the channeled constructs had a large increase of nutrient availability related to the solid counterparts, with an up to 136-fold increase in minimum glucose concentrations and up to 220-fold increase in minimum oxygen concentrations. Under the model assumptions, ECM matrix synthesis increased up to 50% in the constructs with channel, favoring already a small positive impact on the mechanical properties after 72h as a result of improved homogeneity across the tissues (Figure 2).

Fig. 2
Impact of several construct geometrical configurations in free swelling culture conditions on the radial distribution of biphasic mechanical properties after 72h in culture. Reprinted from [6]. This work was expanded to a long-term culture case by simulating an experiment with both a solid and a channeled 2% w/w agarose construct with chondrocytes during 56 days. The model was calibrated with solid construct data on GAG and collagen concentrations, as well as on the compressive Young's modulus and validated by reproducing well the experimental data for the channeled construct [9]. This modeling effort allowed to gain quantitative insights on the spatial heterogeneity of the constructs, showing that the degree of spatial heterogeneity of the Young's modulus the constructs with a central diffusion channel is 23% of the control value, while for permeability the heterogeneity is 27% of the control one, showing a significant improvement for the channeled condition ( Figure 3). The degree of spatial heterogeneity and the insufficient permeability remodeling in the simulated solid constructs affects significantly the mechanical response to compressive strain, with nominal stresses for the simulated heterogeneous TE cartilage 57% lower than for native articular cartilage and pore pressures 53% lower than the native case [9]. Therefore, permeability is the main parameter to be improved to get a more similar mechanical response, c alling for new scaffold material designs and stimulation protocols. The developed complete model was also applied to simulate dynamic loading conditions. A model parametrized with literature parameters was used to simulate the distribution of cell density and ECM in cubic constructs subjected to either compression, shear or bending at 5% of height, 1Hz for 6h continuously [7]. While the simulation time is very short for relevant differences in the mechanical properties to be seen, bending was, under these conditions, the more favourable regime for cell proliferation and the spatial distributions are relevant with the establishment of cart ilage with different structural organizations due to different maximu m principal strain directions ( Figure 4). Current work is related to the model validation for the estimation of cell proliferation, ECM growth and mechanical properties remodelling under several different regimes of cyclic unconfined compressive loading. All the reported works until now have been focused on articular cartilage. However, the complete modeling scheme was also successfully applied to estimate the remodelling of temporomandibu lar joint disc, an area where tissue engineering is still in a very early phase. The application of static hydrostatic pressure for 72h on PEDGA -condylar constructs promoted a very slight improvement of the mechanical properties. This preliminary study pro vided a future basis to estimate the impact of long term static or dynamic hydrostatic pressure on the growth of new temporomandibular joint discs [8] The presented model applications for simulation of growth and remodeling under dynamic loading represent shortterm loading with limited differences in the mechanical properties between conditions. Future work involves the simulation of long-term intermitted compressive loading. On long-term regimes, the collagen content is relevant in terms of fibre organization, therefore the expansion of the constitutive relatio nships to include the anisotropic behaviour of tissue engineered cartilage and its time dynamics is a future goal.

Concluding remarks
This chapter provided an overview of the current state of the computational models used for simulation of growth of tissue engineered cartilage, namely the models based in mixture theory. The several individual contributions for the underlying phenomena, s uch as solute transport and uptake, cell growth, production of extracellular matrix and remodelling of mechanical properties were presented, with a focus on, when applicable, the models that incorporate the impact of mechanical stimulation on these phenomena. Several constitutive relationships to model the mechanical behaviour of tissue engineered cartilage have been proposed and integrated into these models. A fully coupled modeling approach was developed to accommodate all these phenomena in a simultaneou s fashion for a more realistic representation of the biomechanical phenomena and used to estimate the growth and remodelling of mechanical properties under free swelling and mechanical loading of tissues. Future challenges on this area include refinement o f the developed equations through validation with experimental studies with proper measurements of all the underlying variables when possible, study the degree of robustness and/or specificity with different cell-matrix systems, simulation of long term mechanical loading cultures and the incorporation of anisotropic models with reorientation of collagen fibers in the fully coupled formulation.