Continuum Modeling and Simulation in Bone Tissue Engineering

applied sciences Review Continuum Modeling and Simulation in Bone Tissue Engineering Jose A. Sanz-Herrera * and Esther Reina-Romo Escuela Técnica Superior de Ingeniería, Universidad de Sevilla, 41092 Sevilla, Spain *Correspondence: [email protected]; Tel.: +34-95-448-7293 Received: 30 July 2019; Accepted: 23 August 2019; Published: 5 September 2019   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 different stages and problems of interest found in the field 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 classified into: (i) Mechanics modeling and scaffold design, (ii) transport and flow 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 different 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 offers in BTE, and second, the need for interdisciplinary collaboration to further validate numerical models developed in BTE. Keywords: bone tissue engineering; biomaterials; computational mechanobiology; numerical methods in bioengineering 1. 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 Appl. Sci. 2019,9, 3674; doi:10.3390/app9183674 www.mdpi.com/journal/applsci Appl. Sci. 2019,9, 3674 2 of 17 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–10]. On the other hand, polymeric biomaterials are used since they are biodegradable, either natural polymer matrices: Polysaccharides (starch, alginate, chitin/chitosan, hyaluronic acid derivates) or proteins (soy, collagen, fibrin gels, silk). Moreover, synthetic biodegradable polymers for BTE include saturated aliphatic polyesters: Poly(lactic acid) (PLA), poly(glycolic acid) (PGA), poly(lactic-co-glycolide) (PLGA) copolymers, and poly(e-caprolactone); unsaturated linear polyesters: Polypropylene fumarate (PPF); and aliphatic polyesters: Polyhydroxyalkanoates and derived products (PHB, PHBV, P4HB, PHBHHx, PHO). Due to the limited mechanical strength of polymeric biomaterials, a variety of biofibers such as lignocellulosic natural fibers are used as a reinforcement [11]. 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 Appl. Sci. 2019,9, 3674 3 of 17 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. 2. Mechanics Modeling 2.1. 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 – 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 C0as follows, C0=1 VZV Cϕ(Ikh +εϕ(χkh))dV, (1) where εϕ(χkh ) are elementary unit strain solutions with associated displacement field χkh . Ikh is the identity fourth order tensor, Cϕ the stiffness tensor of the base (bulk) material, and Vthe volume of the microscopic cell. On the other hand, the permeability tensor K is obtained following the homogenization theory as, K=1 VZV κi jdV, (2) Appl. Sci. 2019,9, 3674 4 of 17 with κi jthe 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 – 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]. Appl. Sci. 2019, 9, x FOR PEER REVIEW 4 of 16 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–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]. Figure 1. Potential lattice-based scaffold microstructures analyzed for bone tissue engineering (BTE) applications. Reprinted with permission from [37]. 2.2. 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 Figure 1. Potential lattice-based scaffold microstructures analyzed for bone tissue engineering (BTE) applications. Reprinted with permission from [37]. 2.2. 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 Appl. Sci. 2019,9, 3674 5 of 17 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 – 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 – 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 – 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 – 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 Appl. Sci. 2019,9, 3674 6 of 17 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. Appl. Sci. 2019, 9, x FOR PEER REVIEW 6 of 16 Figure 2. Multiscale computer simulation in BTE to predict new bone tissue growth in a unit cell of a predefined scaffold microstructure. (a) CT reconstruction of a rabbit femur, (b) model of the femur, scaffold implantation, loads and boundary conditions, (c) detailed mesh of the model, (d) unit cell microstructure of the scaffold solid domain, and (e) unit cell microstructure of the scaffold fluid domain. Reprinted with permission from [69]. 3. 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–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 Figure 2. Multiscale computer simulation in BTE to predict new bone tissue growth in a unit cell of a predefined scaffold microstructure. ( a ) CT reconstruction of a rabbit femur, ( b ) model of the femur, scaffold implantation, loads and boundary conditions, ( c ) detailed mesh of the model, ( d ) unit cell microstructure of the scaffold solid domain, and ( e ) unit cell microstructure of the scaffold fluid domain. Reprinted with permission from [69]. 3. 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 – 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 Appl. Sci. 2019,9, 3674 7 of 17 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]. Figure 3. Fluid velocity streamlines along poly(L-lactic acid) non regular scaffold computed by means of computational fluid dynamics (CFD) analysis. Reprinted with permission from [84]. 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. 4. 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 Figure 3. Fluid velocity streamlines along poly(l-lactic acid) non regular scaffold computed by means of computational fluid dynamics (CFD) analysis. Reprinted with permission from [84]. 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. Appl. Sci. 2019,9, 3674 8 of 17 4. 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=α∆din Ω+boundary and initial conditions, (3) 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, W=−βdin Ω. (4) 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 ]. Otherresearchers [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–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 Appl. Sci. 2019,9, 3674 9 of 17 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–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). Figure 4. Validation of the simulation of magnesium screw degradation obtained by the mathematical modeling presented by Sanz-Herrera et al. [118] (bottom), compared with the experimental setup by Li et al. [119] (top). Reprinted with permission from [118]. 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. 5. 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–127]. Numerical methods can be particularly useful in this context for the design of the printing strategy and simulation of the biofabrication technique [128–133]. Figure 4. Validation of the simulation of magnesium screw degradation obtained by the mathematical modeling presented by Sanz-Herrera et al. [ 118 ] ( bottom ), compared with the experimental setup by Li et al. [119] (top). Reprinted with permission from [118]. 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. 5. 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 Appl. Sci. 2019,9, 3674 16 of 17 100. Dhote, V.; Vernerey, F.J. Mathematical model of the role of degradation on matrix development in hydrogel scaffold. Biomech. Model. Mechanobiol. 2004,13, 167–183. [CrossRef] 101. Akalp, U.; Bryant, S.J.; Vernerey, F.J. Tuning tissue growth with scaffold degradation in enzyme-sensitive hydrogels: A mathematical model. Soft Matter 2016,12, 7505–7520. [CrossRef] 102. Chen, Y.; Zhou, S.; Li, Q. Mathematical modeling of degradation for bulk-erosive polymers: Applications in tissue engineering scaffolds and drug delivery systems. Acta Biomater. 2011,7, 1140–1149. [CrossRef] 103. Hench, L.L.; Paschall, H.A. Direct chemical bond of bioactive glass ceramic materials to bone and muscle. J. Biomed. Mater. Res. Symp. 1973,4, 25–42. [CrossRef] 104. Wilson, J.; Low, S.B. Bioactive ceramics for periodontal treatment: Comparative studies in the Patus monkey. J. Appl. Biomater. 1992,3, 123–169. [CrossRef] [PubMed] 105. Hench, L.L.; West, J.K. Biological applications of bioactive glasses. Life Chem. Rep. 1996,13, 187–241. 106. Roether, J.A.; Boccaccini, A.R.; Hench, L.L.; Maquet, V.; Gautier, S.; Jerome, R. Development and in vitro characterisation of novel bioresorbable and bioactive composite materials based on polylactide foams and bioglasss for tissue engineering applications. Biomaterials 2002,23, 3871–3878. [CrossRef] 107. Chen, Q.Z.; Thompson, I.D.; Boccaccini, A.R. 45S5 Bioglass derived glass-ceramic scaffolds for bone tissue engineering. Biomaterials 2006,27, 2414–2425. [CrossRef] [PubMed] 108. Sanz-Herrera, J.A.; Boccaccini, A.R. Modelling bioactivity and degradation of bioactive glass based tissue engineering scaffolds. Int. J. Solids Struct. 2011,48, 257–268. [CrossRef] 109. Manhas, V.; Guyot, Y.; Kerckhofs, G.; Chai, Y.C.; Geris, L. Computational modelling of local calcium ions release from calcium phosphate-based scaffolds. Biomech. Model. Mechanobiol. 2017 ,16, 425–438. [CrossRef] [PubMed] 110. Sanz-Herrera, J.A.; Soria, L.; Reina-Romo, E.; Torres, Y.; Boccaccini, A.R. Model of dissolution in the framework of tissue engineering and drug delivery. Biomech. Model. Mechanobiol. 2018 ,17, 1331–1341. [CrossRef] 111. Kraus, T.; Fischerauer, S.F.; Hänzi, A.C.; Uggowitzer, P.J.; Löffler, J.F.; Weinberg, A.M. Magnesium alloys for temporary implants in osteosynthesis: In vivo studies of their degradation and interaction with bone. Acta Biomater. 2012,8, 1230–1238. [CrossRef] 112. Brar, H.S.; Platt, M.O.; Sarntinoranont, M.; Martin, P.I.; Manuel, M.V. Magnesium as a biodegradable and bioabsorbable material for medical implants. JOM 2009,61, 31–34. [CrossRef] 113. Castellani, C.; Lindtner, R.A.; Hausbrandt, P.; Tschegg, E.; Stanzl-Tschegg, S.E.; Zanoni, G.; Beck, S.; Weinberg, A.M. Bone-implant interface strength and osseointegration: Biodegradable magnesium alloy versus standard titanium control. Acta Biomater. 2011,7, 432–440. [CrossRef] 114. Staiger, M.P.; Pietak, A.M.; Huadmai, J.; Dias, G. Magnesium and its alloys as orthopedic biomaterials: A review. Biomaterials 2006,27, 1728–1734. [CrossRef] [PubMed] 115. Grogan, J.A.; Brien, B.J.O.; Leen, S.B.; McHugh, P.E. A corrosion model for bioabsorbable metallic stents. Acta Biomater. 2011,7, 3523–3533. [CrossRef] [PubMed] 116. Grogan,J.A.; Leen,S.B.; McHugh,P.E.Aphysicalcorrosionmodelforbioabsorbablemetalstents. Acta Biomater. 2014,10, 2313–2322. [CrossRef] [PubMed] 117. Bajger, P.; Ashbourn, J.M.A.; Manhas, V.; Guyot, Y.; Lietaert, K.; Geris, L. Mathematical modelling of the degradation behaviour of biodegradable metals. Biomech. Model. Mechanobiol. 2017 ,16, 227–238. [CrossRef] [PubMed] 118. Sanz-Herrera, J.A.; Reina-Romo, E.; Boccaccini, A.R. In silico design of magnesium implants: Macroscopic modeling. J. Mech. Behav. Biomed. Mater. 2018,79, 181–188. [CrossRef] [PubMed] 119. Li, Z.; Sun, S.; Chen, M.; Fahlman, B.D.; Liu, D.; Bi, H. In vitro and in vivo corrosion, mechanical properties and biocompatibility evaluation of MgF 2 -coated MgZn-Zr alloy as cancellous screws. Mater. Sci. Eng. C 2017,75, 1268–1280. [CrossRef] [PubMed] 120. Frenning, G.; Brohede, U.; Stromme, M. Finite element analysis of the release of slowly dissolving drugs from cylindrical matrix systems. J. Control. Release 2005,107, 320–329. [CrossRef] 121. Peppas, N.A.; Narasimhan, B. Mathematical models in drug delivery: How modeling has shaped the way we design new drug delivery systems. J. Control. Release 2014,190, 75–81. [CrossRef] 122. Mozafari, M.; Sefat, F.; Atala, A. Handbook of Tissue Engineering Scaffolds; Elsevier: Amsterdam, The Netherlands, 2019; Volume 1, ISBN 9780081025642. Appl. Sci. 2019,9, 3674 17 of 17 123. Visconti, R.P.; Kasyanov, V.; Gentile, C.; Zhang, J.; Markwald, R.R.; Mironov, V. Towards organ printing: Engineering an intra-organ branched vascular tree. Expert. Opin. Biol. Ther. 2010,10, 409–420. [CrossRef] 124. Derby, B. Printing and prototyping of tissues and scaffolds. Science 2012 ,338, 921–992. [CrossRef] [PubMed] 125. Murphy, S.V.; Atala, A. 3D bioprinting of tissues and organs. Nat. Biotechnol. 2014 ,32, 773–785. [CrossRef] [PubMed] 126. Dababneh, A.B.; Ozbolat, I.T. Bioprinting technology: A current state-of-the-art review. J. Manuf. Sci. Eng. 2014,136, 061016. [CrossRef] 127. Mandrycky, C.; Wang, Z.; Kim, K.; Kim, D.H. 3D bioprinting for engineering complex tissues. Biotechnol. Adv. 2016,34, 422–434. [CrossRef] [PubMed] 128. Li, M.; Tian, X.; Schreyer, D.J.; Chen, X. Effect of needle geometry on flow rate and cell damage in the dispensing-based biofabrication process. Biotechnol. Prog. 2011,27, 1777–1784. [CrossRef] 129. Guyot, Y.; Papantoniou, I.; Chai, Y.C.; Van Bael, S.; Schrooten, J.; Geris, L. A computational model for cell/ECM growth on 3D surfaces using the level set method: A bone tissue engineering case study. Biomech. Model. Mechanobiol. 2014,13, 1361–1371. [CrossRef] [PubMed] 130. Guyot, Y.; Luyten, F.P.; Schrooten, J.; Papantoniou, I.; Geris, L. A three dimensional computational fluid dynamics model of shear stress distribution during neotissue growth in a perfusion bioreactor. Biotechnol. Bioeng. 2015,112, 2591–2600. [CrossRef] 131. Guyot, Y.; Papantoniou, I.; Luyten, F.P.; Geris, L. Coupling curvature dependent and shear stress-stimulated neotissue growth in dynamic bioreactor cultures: A 3D computational model of a complete scaffold. Biomech. Model. Mechanobiol. 2016,15, 169–180. [CrossRef] 132. Billiet, T.; Gevaert, E.; De Schryver, T.; Cornelissen, M.; Dubruel, P. The 3D printing of gelatin methacrylamide cell-laden tissue-engineered constructs with high cell viability. Biomaterials 2014,35, 49–62. [CrossRef] 133. Reina-Romo, E.; Papantoniou, I.; Bloemen, V.; Geris, L. Computational design of tissue engineering scaffolds. In Handbook of Tissue Engineering Scaffolds; Elsevier: Amsterdam, The Netherlands, 2019. © 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).