Continuum Modeling and Simulation in Bone Tissue Engineering

Featured Application: Bone tissue engineering (BTE) can be investigated by means of mathematical modeling and numerical tools complementary to the experimental methods. In particular, the design process of a bone tissue engineering product or protocol can be enhanced by computer simulation as evidenced in a number of papers in the last decades. In this work, we review the most relevant contributions of continuum models and simulations applied to di ﬀ erent stages and problems of interest found in the ﬁeld of BTE. Abstract: Bone tissue engineering is currently a mature methodology from a research perspective. Moreover, modeling and simulation of involved processes and phenomena in BTE have been proved in a number of papers to be an excellent assessment tool in the stages of design and proof of concept through in-vivo or in-vitro experimentation. In this paper, a review of the most relevant contributions in modeling and simulation, in silico, in BTE applications is conducted. The most popular in silico simulations in BTE are classiﬁed into: (i) Mechanics modeling and sca ﬀ old design, (ii) transport and ﬂow modeling, and (iii) modeling of physical phenomena. The paper is restricted to the review of the numerical implementation and simulation of continuum theories applied to di ﬀ erent processes in BTE, such that molecular dynamics or discrete approaches are out of the scope of the paper. Two main conclusions are drawn at the end of the paper: First, the great potential and advantages that in silico simulation o ﬀ ers in BTE, and second, the need for interdisciplinary collaboration to further validate numerical models developed in BTE.


Introduction
Bone tissue engineering (BTE) aims to persuade bone tissue to regenerate and/or heal under diverse circumstances. These circumstances include the repair of long bone defects where the body has limited regeneration capability, treatment of bone diseases such as osteoporosis, or to accelerate the process of bone fracture healing [1]. The bone healing process can be summarized in the following steps [2]: (i) A hematoma is formed from injury in the periosteum; (ii) osteocytes near the fracture site die due to blood disruption following a demand for repair; (iii) macrophages and fibroblasts are recruited to the site to remove tissue debris and to express extracellular matrix, then growth factors and cytokines released by these inflammatory cells, mesenchymal stem cells are recruited from the bone marrow and periosteum, which then proliferate and differentiate into progenitor cells; (iv) those osteoprogenitor cells differentiate into osteoblasts and form osteoid which is rapidly calcified into bone. Finally, (v) the uncalcified material is resorbed and new bone is deposited. The woven bone is then remodeled into lamellar bone and the process is completed by the return of normal bone marrow within trabecular regions, while in repairing cortical bone the spaces between trabeculae are gradually filled in with successive layers of bone thus forming new Haversian canals.
In order to mimic the process of natural bone healing described above, BTE methodology involves the use of porous biomaterials as a structure support, i.e., scaffolds, which serve to temporally cover the bone defect, as well as providing room for bone cells to invade, proliferate, and develop their specific functions to segregate new bone matrix. Ideally, the structural support (scaffold) should degrade over time resulting in the formation of a new bony structure [1].
The base scaffold biomaterial should mimic the natural bone matrix; therefore, ceramics and polymers are extensively used in BTE applications. On the one hand, bioceramics including hydroxyapatite (HA) and calcium phosphate are used. They show excellent ability to bond to bone, good biocompatibility behavior, and reasonably good mechanical properties [3][4][5][6][7][8][9][10].
Bioactive glasses are a class of inorganic biomaterials discovered by Hench in 1969 [12]. Bioactive glass materials have the ability to bond new bone tissue enhanced by a reaction that takes place once the biomaterial is inserted in the body environment. A number of commercial scaffolds have been developed such as Bioglass ® in bone defects, PerioGlas TM for periodontal disease, and NovaBone TM as a bone filler [13]. Some drawbacks of bioactive glasses are their low fracture toughness and mechanical strength, especially in a porous form. These drawbacks are enhanced by a HA particulate reinforcement [14].
In order to promote new bone tissue regeneration, following the natural process of bone healing, the BTE methodology usually combines a seeding strategy along the biomaterial surface of the scaffold prior to implantation. Usually bone marrow stromal cells (BMSC), mesenchymal stem cells (MSC), or preosteoblasts are combined with bone morphogenetic proteins (BMP) or growth factors such as transforming growth factor-β (TGF-β), within porous scaffolds with the use of a bioreactor. In particular, BMSC-seeded ceramics were implanted in a large tibial defect in ewes [15]. Mechanical stability at the defect site was obtained either by an internal plate or by external fixation. Results show around 10% of bone volume to total volume regeneration. Mastrogiacomo et al. [16] also presented the same tissue engineering strategy for segmental tibial defect in sheep, resulting in a progressive scaffold resorption coincident with new bone deposition [16]. A comparison of bone regeneration in rabbits among BMSC-harvested poly(lactide-co-glycolide) scaffolds, non-harvested ones, and without scaffolds was presented in [17]. The defect consisted of a unilateral femoral osteotomy gap created surgically under general anesthetic and stabilized by a mandibular reconstruction plate fixated with three screws on either side of the osteotomy site. The obtained results showed that there were no significant differences between the use of the scaffold alone and the cell-seeded scaffold although faster in this case. On the contrary, the non-scaffold strategy presented less regenerated tissue. Using biodegradable scaffolds as well, Holy et al. [18] performed a trabecular-like, three-dimensional structure to repair bone. BMCS-preseeded scaffolds were implanted in a non-healing rabbit segmental bone defect achieving bony union within eight weeks. Biodegradable scaffolds medical-grade polycaprolactone-calcium phosphate (mPCL-CaP) seeded with BMSC were implanted under the skin of nude rats. It presented neo cortical and well-vascularized cancellous bone up to 40% of bone volume [19]. Savarino et al. [20] presented a study involving the use of poly-e-caprolactone scaffolds (PCL) loaded with BMSC and BMP-4. PCL without cells showed scarce bone formation and scaffold resorption, whereas PCL seeded with BMSC stimulated new tissue formation. In conclusion, the combination of BMSC with BMP-4 strongly favored osteoinductivity of cellular constructs.
Bone tissue undergoes a piezoelectric effect such that bone cells are stimulated to remodel by means of an electrical filed, induced by a strain field, which recruits and aggregates macromolecules and ions in the extracellular matrix [21]. This feature has been exploited in bone tissue regeneration using piezoelectric scaffolds which mimic the natural collagen matrix, such as piezopolymers. The reader is referred to [22] for a review of piezoelectric scaffolds.
Mathematical modeling and computer simulation have been demonstrated to be a powerful tool both in the design and evaluation stages of scaffolds in BTE applications. In the design phase, the mechanical properties (elasticity modulus, strength, toughness) and fluidic properties (diffusivity and permeability) are essential for the overall success of the scaffold. Moreover, the degradation properties and characteristics are also of interest for a certain BTE experiment. The overall behavior during implantation of BTE processes has also been simulated through the incorporation of biomechanical and mechanobiological theories. In the next sections, the state-of-the-art of continuum approaches and simulation of the referred BTE phenomena are reviewed. The new challenges, perspectives, and some conclusions are drawn at the end of the paper.

Constitutive Behavior Modeling and Scaffold Design
According to the evidences, the optimal scaffold design, i.e., microarchitecture, should include high porosity, proper pore structure interconnection, and enough specific surface to attach cells to segregate new matrix and proliferate [1]. In particular, porosity values in the range 60-80% are recommended for load-bearing BTE applications [1,11]. These features can be controlled via a computer aided design (CAD) technique of the scaffold in the design step. On the other hand, permeability is an important parameter in tissue engineering applications, linked to pore structure and porosity, and related to the flow of nutrients and waste removal which are essential processes during cell activity. Moreover, the overall constitutive mechanical behavior is a design variable of the scaffold microstructure for load bearing in BTE, as well as the distribution of the mechanical stimuli to activate bone cell functions [23].
Both permeability and overall mechanical properties are macroscopic quantities which describe the fluid and solid mechanical behavior of the scaffold. They can be tailored by getting control over some scaffold parameters microscopically, e.g., mean pore size, porosity, and virgin biomaterial mechanical properties. Hollister and co-workers have focused on the design of bone scaffolds as an optimization problem to get a microstructure as similar as possible (mechanical properties, porosity, pore size, etc.) to that of the implanted region. For this purpose, the homogenization theory was extensively applied during design [24][25][26]. In the same context, the asymptotic homogenization theory was applied to validate the experimental Darcian permeability and mechanical properties of a specific scaffold, Sponceram ® , available for BTE applications [27]. On the one hand, the homogenization theory computes the overall macroscopic stiffness tensor C 0 as follows, where ε ϕ (χ kh ) are elementary unit strain solutions with associated displacement field χ kh . I kh is the identity fourth order tensor, C ϕ the stiffness tensor of the base (bulk) material, and V the volume of the microscopic cell.
On the other hand, the permeability tensor K is obtained following the homogenization theory as, with κ i j the characteristic fluid velocities associated with unit pressure gradients. The homogenization theory has been further developed in connection with topology optimization in the scaffold design [28]. Coelho et al. [29] used again the homogenization theory and topological optimization connected to the fabrication and testing of the optimal scaffold. Moreover, the mechanical properties of a multiphasic scaffold composed of poly(propylene fumarate) reinforced with silicon particles were obtained in terms of elasticity and shear moduli by means of theoretical developments based on the Eshelby theory [30]. The overall mechanical properties of porous and multiphasic scaffolds have been also obtained by means of theoretical (analytical) methods [31][32][33]. A review of these methods can be seen in [34].
As seen in the cited papers above, the design of the scaffold microstructure can be posed as an optimization problem to achieve the optimal overall constitutive behavior. The same optimization design procedure was followed to achieve a uniform shear stress as an objective function [35]. This criterion was established to get an optimal erosion along the microstructure. In the same context, Boccaccio et al. [36] determined the optimal loading for a number of 3D printed scaffold microarchitectures to enhance new bone tissue growth. This study may be used in patient-specific BTE applications.
The design of the scaffold microstructure has been recently linked to advanced manufacturing techniques mainly based on 3D printing. Egan et al. [37] analyzed lattice-based scaffold microstructural parameters of different families ( Figure 1). Results concluded that each family may find an application depending on its specific characteristics. The mechanical performance of 3D printed scaffolds was evaluated as well with a focus on the distribution of the mechanical stimuli along the scaffold microstructure [38]. The homogenization theory has been further developed in connection with topology optimization in the scaffold design [28]. Coelho et al. [29] used again the homogenization theory and topological optimization connected to the fabrication and testing of the optimal scaffold. Moreover, the mechanical properties of a multiphasic scaffold composed of poly(propylene fumarate) reinforced with silicon particles were obtained in terms of elasticity and shear moduli by means of theoretical developments based on the Eshelby theory [30]. The overall mechanical properties of porous and multiphasic scaffolds have been also obtained by means of theoretical (analytical) methods [31][32][33]. A review of these methods can be seen in [34].
As seen in the cited papers above, the design of the scaffold microstructure can be posed as an optimization problem to achieve the optimal overall constitutive behavior. The same optimization design procedure was followed to achieve a uniform shear stress as an objective function [35]. This criterion was established to get an optimal erosion along the microstructure. In the same context, Boccaccio et al. [36] determined the optimal loading for a number of 3D printed scaffold microarchitectures to enhance new bone tissue growth. This study may be used in patient-specific BTE applications.
The design of the scaffold microstructure has been recently linked to advanced manufacturing techniques mainly based on 3D printing. Egan et al. [37] analyzed lattice-based scaffold microstructural parameters of different families ( Figure 1). Results concluded that each family may find an application depending on its specific characteristics. The mechanical performance of 3D printed scaffolds was evaluated as well with a focus on the distribution of the mechanical stimuli along the scaffold microstructure [38].

Simulation of Applications of Interest
In silico modeling in BTE is becoming a predictive tool for bone regeneration simulation within the scaffold reducing the number of in vivo/ in vitro experiments that need to be performed.
Bone tissue regeneration in vivo using scaffolds inherently attends to two well-differentiated spatial and temporal scales. One is the tissue level or macroscopic scale, with the other one being named as pore level or micro/mesoscopic scale. Therefore, the mathematical and in silico modeling of tissue growth within scaffolds has been usually restricted to one of these scales.
At the macroscopic scale, finite element analyses (FEA) has been used to analyze the solid mechanical behavior of the scaffold. Two different modeling approaches can be found in the

Simulation of Applications of Interest
In silico modeling in BTE is becoming a predictive tool for bone regeneration simulation within the scaffold reducing the number of in vivo/in vitro experiments that need to be performed.
Bone tissue regeneration in vivo using scaffolds inherently attends to two well-differentiated spatial and temporal scales. One is the tissue level or macroscopic scale, with the other one being named as pore level or micro/mesoscopic scale. Therefore, the mathematical and in silico modeling of tissue growth within scaffolds has been usually restricted to one of these scales.
At the macroscopic scale, finite element analyses (FEA) has been used to analyze the solid mechanical behavior of the scaffold. Two different modeling approaches can be found in the literature at this scale: FEA models, that simulate the fluid flow with poroelastic materials (in which the flow is defined according to Darcy's law) and computational fluid dynamics (CFD) models, which solve the governing equations of fluid flow dynamics.
On the one hand, FEA have tried to investigate the effect of scaffold architecture on the neotissue formed within the scaffold. The mathematical model of cell/tissue differentiation proposed by Prendergast et al. [39] has been widely used in bone tissue engineering applications [40][41][42][43][44][45][46][47]. In these studies, cell proliferation is modeled as a diffusion problem and macroscopic variables, such as the shear strain and fluid flow, are considered as stimuli for cell differentiation. Osteochondral defect repair using scaffolds has been simulated in [40]. Khayyeri et al. [42] predicted tissue differentiation with FEA by means of a lattice modeling approach in a bone chamber and compared the numerical results with in vivo bone formation. Olivares et al. [41] used a CAD-based model to analyze the effects of scaffold pore morphology on cell differentiation stimulus values. These models are either based on computer aided design (CAD) [42] or on micro-computed tomography (µCT) [43,44,47]. However, CAD models are not able to capture the strain distributions seen in the µCT [48].
On the other hand, CFD characterizes the hydrodynamic field imposed on cells within different types of scaffolds [48][49][50][51][52][53]. They allow numerical prediction of the shear stresses induced at the internal walls in relation with the scaffold geometry and, thus, they qualitatively relate the flow of culture medium with the macroscopic 3D shear stresses [54,55]. It is well known that flow-induced shear stress produces a significant stimulatory effect and plays a critical role in the physiological responses [56][57][58][59]. More recent studies have predicted the new tissue formed and have related the wall shear stresses (WSS) with cellular processes such as proliferation, differentiation, or mineralization of the extracellular matrix by means of CFD analyses. In particular, Sonnaert et al. [60] investigated the influence of fluid-flow-induced shear stress on the proliferation, differentiation, and matrix deposition of human periosteal-derived cells in 3D Ti6Al4V scaffolds. Nava et al. [61] numerically described the growth of a neotissue in a perfusion bioreactor where the growth was coupled to oxygen concentration and shear stress. Zhao et al. [62] calculated the optimal fluid flow rate to be applied to the bioreactor for BTE experiments by means of CFD for different scaffold pore shapes, pore diameters, and porosities. They employed a combination of CFD with mechano-regulation theories to optimize the external flow rate for a perfusion bioreactor. Such flow rate would maximize the scaffold surface fraction, whose WSS was in the range required for mineralization.
In the microscopic scale, Byrne et al. [63] have presented a mechanobiological model applied to a periodic unit cell of the scaffold microstructure. In this model, the influence of several factors such as permeability, mechanical properties, and others have been analyzed in tissue regeneration. In the same scale, bone tissue regeneration within a scaffold unit cell has been numerically reproduced [64] using concepts and hypothesis previously established for a remodeling theory [65]. At a nanoscale level, the mechanisms of cell adhesion to the walls has also been studied [66,67].
Finally, the use of multiscale methods and homogenization has allowed the analysis of different phenomena of bone regeneration within scaffolds [68][69][70][71][72][73]. Sanz-Herrera et al. [68,69] proposed for the first time a coupled micro-macro mathematical approach for bone tissue regeneration in tissue engineering applications, see Figure 2. In these works, at the tissue level, the macroscopic mechanical, diffusive, and flow properties are derived by means of the asymptotic homogenization theory. At the microscopic scale, bone tissue regeneration at the scaffold microsurface is simulated using a bone growth model based on a bone remodeling theory [23], whereas scaffold degradation is implemented following a previous model [64]. This multiscale model was later used to analyze the effect of scaffold microarchitecture in new bone tissue growth [72]. Nguyen et al. [74] performed a multiscale approach based on a CFD analysis on two scales to predict cell growth inside a given macroporous scaffold put in a perfusion bioreactor. Their objective was to determine the optimal flow rate in order to enhance cell proliferation and to improve an upcoming bone reconstruction.

Transport and Flow Modeling
Fluid circulation, transport of nutrients, and waste removal are essential processes that must take place inside the scaffold microstructure in a successful BTE product. In fact, vascularization, i.e., invasion of blood vessel and formation of new vasculature, is one of the main drawbacks in tissue engineering [75]. Moreover, fluid percolation and circulation is important as well in vitro in BTE assisted by bioreactors. Therefore, flow mechanics and diffusion mechanisms have been theoretically modeled and simulated during design of BTE scaffolds.
Permeability is an overall parameter that (macroscopically) condenses the information about how a fluid penetrates in a porous medium. It is a clear design parameter in BTE scaffolds. Scaffold permeability has been measured experimentally in several papers using pumped water [76,77], with a gravity-induced pressure using water [78][79][80] or using compressed air [81], to cite a few. Indeed, permeability was numerically validated using the homogenization theory versus an experimental gravity-induced pressure setup, for a specific commercial scaffold for BTE applications [27]. In the same context, Truscello et al. [82] validated the permeability values of Ti6Al4V scaffolds using CFD simulations. An interesting study [83] obtained the permeability of simulated cancellous bone structure by means of CFD analyses over a unit cell of the scaffold. Results concluded which wasthe best-suited scaffold microarchitecture for BTE in terms of blood circulation inside the microstructure.
The study presented in [84], see Figure 3, shows a CFD analysis over a complex non regular microstructure of a poly(L-lactic acid) scaffold. The microstructural geometry was obtained by

Transport and Flow Modeling
Fluid circulation, transport of nutrients, and waste removal are essential processes that must take place inside the scaffold microstructure in a successful BTE product. In fact, vascularization, i.e., invasion of blood vessel and formation of new vasculature, is one of the main drawbacks in tissue engineering [75]. Moreover, fluid percolation and circulation is important as well in vitro in BTE assisted by bioreactors. Therefore, flow mechanics and diffusion mechanisms have been theoretically modeled and simulated during design of BTE scaffolds.
Permeability is an overall parameter that (macroscopically) condenses the information about how a fluid penetrates in a porous medium. It is a clear design parameter in BTE scaffolds. Scaffold permeability has been measured experimentally in several papers using pumped water [76,77], with a gravity-induced pressure using water [78][79][80] or using compressed air [81], to cite a few. Indeed, permeability was numerically validated using the homogenization theory versus an experimental gravity-induced pressure setup, for a specific commercial scaffold for BTE applications [27]. In the same context, Truscello et al. [82] validated the permeability values of Ti6Al4V scaffolds using CFD simulations. An interesting study [83] obtained the permeability of simulated cancellous bone structure by means of CFD analyses over a unit cell of the scaffold. Results concluded which wasthe best-suited scaffold microarchitecture for BTE in terms of blood circulation inside the microstructure.
The study presented in [84], see Figure 3, shows a CFD analysis over a complex non regular microstructure of a poly(l-lactic acid) scaffold. The microstructural geometry was obtained by microCT and results analyzed both the permeability parameter as well as the wall shear stress, which is an important mechanobiological output as a mechanical stimulus [39], as seen before. The experimental results validated the simulations in the referred work. Similar studies were conducted for irregular pore geometries and regular ones [41,46]. Homogeneous fluid flow distribution was found for irregular microstructures, in contrast to regular ones, due to irregular interconnection between the pores [85].
Appl. Sci. 2019, 9, x FOR PEER REVIEW 7 of 16 microCT and results analyzed both the permeability parameter as well as the wall shear stress, which is an important mechanobiological output as a mechanical stimulus [39], as seen before. The experimental results validated the simulations in the referred work. Similar studies were conducted for irregular pore geometries and regular ones [41,46]. Homogeneous fluid flow distribution was found for irregular microstructures, in contrast to regular ones, due to irregular interconnection between the pores [85]. On the other hand, diffusion and transport of nutrients have been simulated in the design of scaffold microstructure. The macroscopic diffusivity of porous scaffolds, as an analogy to Darcian permeability, was analyzed as a function of the pore microarchitecture [68]. These results were used in a full multiscale model accounting for microstructural growth and resorption, and subsequent macroscopic evolution. Li et al. [86] showed finite element simulations of the diffusion and evolution of oxygen concentration along the interior of regular scaffold microstructures.
The evolution of tissue growth coupled with nutrient supply and also including mechanistic effects has been mathematically modeled mostly for in vitro (bioreactor) conditions. In particular, Lemon et al. [87] proposed a continuum and multiphasic model including scaffold, cells, and water (medium) as the main phases. The model was applied to the study of the mobility and aggregation of a population of cells seeded into an artificial polymeric scaffold. This model was extended in [88]. A similar approach, which includes coupling between tissue growth and nutrients consumption, was proposed in [89]. The model was numerically elaborated, implemented, and solved in a hollow-fiber membrane bioreactor as an application, validating existing results available in the literature.

Modeling of Physical Phenomena
In this section, we review the modeling and simulation of the most relevant physical phenomena that take place in biomaterials in BTE applications, such as biodegradation. Most BTE products exploit the ability of biomaterials to dissolve over time, finally resulting in the removal of the scaffold implant. This class of biomaterials are synthetic polymers: Lactide polymer On the other hand, diffusion and transport of nutrients have been simulated in the design of scaffold microstructure. The macroscopic diffusivity of porous scaffolds, as an analogy to Darcian permeability, was analyzed as a function of the pore microarchitecture [68]. These results were used in a full multiscale model accounting for microstructural growth and resorption, and subsequent macroscopic evolution. Li et al. [86] showed finite element simulations of the diffusion and evolution of oxygen concentration along the interior of regular scaffold microstructures.
The evolution of tissue growth coupled with nutrient supply and also including mechanistic effects has been mathematically modeled mostly for in vitro (bioreactor) conditions. In particular, Lemon et al. [87] proposed a continuum and multiphasic model including scaffold, cells, and water (medium) as the main phases. The model was applied to the study of the mobility and aggregation of a population of cells seeded into an artificial polymeric scaffold. This model was extended in [88]. A similar approach, which includes coupling between tissue growth and nutrients consumption, was proposed in [89]. The model was numerically elaborated, implemented, and solved in a hollow-fiber membrane bioreactor as an application, validating existing results available in the literature.

Modeling of Physical Phenomena
In this section, we review the modeling and simulation of the most relevant physical phenomena that take place in biomaterials in BTE applications, such as biodegradation. Most BTE products exploit the ability of biomaterials to dissolve over time, finally resulting in the removal of the scaffold implant. This class of biomaterials are synthetic polymers: Lactide polymer (trimethylene carbonate d,l-lactide (TMCDLLA)) [90], poly-e-caprolactone (PCL) [91,92], polylactic acid (PLA) [93], and polyglycolic acid (PGA) [94,95], to cite a few. Even though the degradation products are naturally evacuated from the human body, some secondary problems have been reported in the reaction of these residuals with the living tissues [96].
Modeling and simulation of biodegradation of polymeric biomaterials were mainly developed after the theory proposed by Gopferich [97]. This modeling was the basis for the development of a number of applications in BTE using polymeric scaffolds and it considers water diffusion within the bulk polymeric biomaterial according to a Fickean law as follows, . d = α∆d in Ω + boundary and initial conditions, where d is the water (aqueous) concentration, α the diffusion coefficient, and ∆ the Laplacian operator. The dot on d denotes time derivative and Ω the biomaterial domain. Water concentration, W, within the polymer is then related to the change rate of the molecular weight of the biomaterial due to hydrolysis and assumed to depend on the local water content, where β is a constant property of the material. The model introduced above was applied to the study of the degradation of a unit cell of the scaffold microstructure [66] and coupled with a bone growth multiscale model [68,69,73]. Other researchers [98,99] presented a similar but extended model to analyze degradation of biopolymers, which takes into consideration crystallization. Further studies [100,101] introduced a mathematical continuum formulation based on the mixture theory to model degradation, transport of molecules across the extracellular matrix, and swelling in a hydrogel scaffold. The model was useful to investigate hydrolytic and enzymatic degradation. Moreover, a similar model for the simulation of the hydrolysis phenomena in polymeric biomaterials was presented in [102]. In this approach, the hydrolysis reaction was modeled by a fundamental stochastic process and an additional autocatalytic effect. Results were shown over different polymeric matrices providing a good agreement with experimental data in the literature.
Inorganic bioactive materials, such as bioglasses, are widely used in BTE due to their ability to react with the body fluid and dissolve, finally resulting in the formation of a hydroxyapatite surface layer. This layer can form stable bonds with the adjacent living tissue, which is specially well suited to bone implants and BTE [11,[103][104][105][106][107].
Dissolution of bioactive glasses was macroscopically modeled using a physical-chemical model which turned into a reaction-diffusion continuum model [108]. Dissolution was represented using the Voxel-FEM approach. Additionally, a reaction-diffusion modeling approach was proposed for the simulation of degradation of calcium phosphate scaffolds [109]. Recently, a generic mathematical framework for the simulation and design of dissolution of biomaterials for tissue engineering and drug delivery applications has been introduced [110]. The model was experimentally validated by means of a straightforward ad hoc setup that considered the dissolution of bicarbonate pellets.
Even though magnesium implants have been used since several decades ago, they have gained an increasing interest in the last years. The main motivations are the great advantages they offer versus traditional metallic, bioceramic, or polymeric implants. On one hand, the absence of a second surgery for removal due to their biodegradability characteristic [111]. Moreover, magnesium implants show similar mechanical properties to bone tissue, avoiding bone resorption due to stress shielding in the neighborhood of the implant and hence minimizing osseointegration problems. Finally, magnesium implants can be used as a load-bearing scaffold being then an ideal candidate for BTE applications [111][112][113]. As a main drawback, magnesium biomaterials release hydrogen gas as a consequence of biodegradation which may prove dangerous under uncontrolled fast dissolution reaction rates [114]. Modeling and simulation of magnesium biodegradation is therefore identified as a useful assessment tool in the design phase of magnesium implants and BTE applications both to understand and control the dissolution process of such biomaterials. In this context, Grogan et al. [115,116] presented physically based continuum models, based on reaction-diffusion equations, available to analyze the corrosion of magnesium metal stents. A similar mathematical modeling was shown in [117] but using the level-set strategy to simulate the dissolution of the biomaterial. Sanz-Herrera et al. [118] developed a continuum model for the biodegradation of magnesium accounting for the dynamics and evolution of secondary species, such as pH or corrosion products. The model was qualitatively validated with previous results found in the literature (Figure 4).
Appl. Sci. 2019, 9, x FOR PEER REVIEW 9 of 16 similar mathematical modeling was shown in [117] but using the level-set strategy to simulate the dissolution of the biomaterial. Sanz-Herrera et al. [118] developed a continuum model for the biodegradation of magnesium accounting for the dynamics and evolution of secondary species, such as pH or corrosion products. The model was qualitatively validated with previous results found in the literature (Figure 4). The general continuum modeling framework proposed above based on reaction-diffusion equations, numerical implementation, and simulation has been applied to drug delivery systems [120]. Additional mathematical models have been proposed in this field in the last decades. Since this is not the focus of this paper, the reader is addressed to [121] for a review of models in drug delivery.

Perspectives
Continuum modeling and simulation in BTE has been proved to be a useful tool at several stages of the methodology according to the revised literature. On one hand, the design of the scaffold in terms of its desired characteristics, such as porosity, pore connectivity, permeability, or overall mechanical behavior, according to a specific application, has been a niche of research using numerical methods. The literature is nowadays vast regarding computer-aided scaffold design, and, traditionally, researchers involved in this part have not always been connected with the fabrication process, which has been performed in collaboration with biomaterial scientists. Currently, 3D printing has brought both fields together, and now computer design of scaffolds is usually The general continuum modeling framework proposed above based on reaction-diffusion equations, numerical implementation, and simulation has been applied to drug delivery systems [120]. Additional mathematical models have been proposed in this field in the last decades. Since this is not the focus of this paper, the reader is addressed to [121] for a review of models in drug delivery.

Perspectives
Continuum modeling and simulation in BTE has been proved to be a useful tool at several stages of the methodology according to the revised literature. On one hand, the design of the scaffold in terms of its desired characteristics, such as porosity, pore connectivity, permeability, or overall mechanical behavior, according to a specific application, has been a niche of research using numerical methods. The literature is nowadays vast regarding computer-aided scaffold design, and, traditionally, researchers involved in this part have not always been connected with the fabrication process, which has been performed in collaboration with biomaterial scientists. Currently, 3D printing has brought both fields together, and now computer design of scaffolds is usually complemented by 3D printing of prototypes and mechanical testing of the specimens in the lab [122]. In addition, bioprinting techniques, i.e., 3D printing technology using biological materials, are being explored as an alternative to traditional BTE [123][124][125][126][127]. Numerical methods can be particularly useful in this context for the design of the printing strategy and simulation of the biofabrication technique [128][129][130][131][132][133].
Other physical phenomena that are relevant to analyze in the design phase of a BTE product have not been treated and investigated in detail by computer simulation. These phenomena include biomaterial dissolution and biodegradation since a limited number of mathematical models have been presented in the literature, despite the clear advantages that introduce simulation-based design of degradable materials as concluded in the referred examples. This is a clear line for future research in the modeling of BTE.
On the other hand, the evolution of the scaffold in terms of tissue growth, tissue differentiation, and tissue-scaffold interaction has been extensively modeled and simulated in BTE applications in the last decades. Some elaborated models including multiscale, multiphasic, and multiphysics elements have been presented in several examples of interest, showing a great potential. First, the models predict in silico the performance of a certain in vivo applications by means of a virtual assay conducted in the computer. The simulation can be patient-specific accounting for the specific characteristics in hands. Second, in silico simulations reduce time and cost of the assays, as well as reducing animal experimentation with its subsequent ethical implications. However, the main drawback of the computer simulation of clinical BTE applications is the fact that the models are based on empirical laws of tissue growth and evolution. These theories have not been sufficiently validated although some of them were proposed several decades ago. It is therefore necessary to advance the validation of the models in connection with a multidisciplinary team of clinicians, veterinarians and biologists, by no means a trivial task. Model validation is then identified as another research line in BTE collateral to in silico BTE. Moreover, multiscale and multiphasic analyses have been limited by the computational cost of this kind of modeling. In particular, multiphasic models traditionally consider both the solid and fluid domains of the scaffold in a simplified way, either using diffusive or poroelastic approaches. A coupled CFD (fluid phase) and mechanical (solid phase) approach, including solid-fluid interaction, can be useful to account in a detailed fashion both the mechanobiological stimulus, as well as degradation phenomena that take place in the biomaterial.
Finally, even a model available for BTE simulation that has been validated at some conditions, may not be predictive in a different scenario. This fact is due to the highly phenomenological nature of this kind of continuum model. Usually, models for BTE include an average of 10-15 model parameters to be calibrated by experimental setups. It is therefore needed to establish a magnified and more fundamental observation scale in the continuum models, even lower than the pore scaffold scale; which allows to consider in a detailed fashion cell-biomaterial sensing and interaction, as well as specific functions of the cell such as proliferation, migration, or apoptosis. These models, posed at the cell scale, may violate the continuum assumption and other discrete models, such as agent-based models, including lattice or particle methods can be of application at this scale. Nonetheless, the discrete models defined at the cell scale may be useful at meso-or macroscopic continuum scale models by fitting complex correlations at this scale a priori or feeding with information at higher scales using a multiscale coupling approach.

Conclusions
Through this review, the usefulness and potential of continuum models and computer simulation, both in the in silico design of scaffolds for BTE and in the prediction of the evolution of bone tissue growth in BTE, have been demonstrated. Even though some issues regarding in silico design of scaffolds can be explored in more detail, these simulations are usually calibrated and contrasted with experimental results. On the other hand, the predictive capability of in silico models of applications of interest in BTE needs to be enhanced. The empirical mechanobiological rules underlying those models need further validation via interdisciplinary collaboration among the different involved researchers. In the improvement of the predictive capability of continuum approaches in BTE, it may help to pose the problem in a higher and magnified observation scale, and hence less phenomenological scale, using discrete numerical approaches such as agent-based modeling or molecular dynamics simulations. All these efforts may contribute to make BTE a clinical viable reality in the next years.
Author Contributions: J.A.S.-H. and E.R.-R. analyzed the state of the art in the field, discussed the previous work and perspectives, and wrote the review paper.
Funding: This research was funded by the Ministerio de Economía y Competitividad del Gobierno España, grant number DPI2017-82501-P.

Conflicts of Interest:
The authors declare no conflict of interest.