Geometry Design Optimization of Functionally Graded Scaffolds for Bone Tissue Engineering: A Mechanobiological Approach

Functionally Graded Scaffolds (FGSs) are porous biomaterials where porosity changes in space with a specific gradient. In spite of their wide use in bone tissue engineering, possible models that relate the scaffold gradient to the mechanical and biological requirements for the regeneration of the bony tissue are currently missing. In this study we attempt to bridge the gap by developing a mechanobiology-based optimization algorithm aimed to determine the optimal graded porosity distribution in FGSs. The algorithm combines the parametric finite element model of a FGS, a computational mechano-regulation model and a numerical optimization routine. For assigned boundary and loading conditions, the algorithm builds iteratively different scaffold geometry configurations with different porosity distributions until the best microstructure geometry is reached, i.e. the geometry that allows the amount of bone formation to be maximized. We tested different porosity distribution laws, loading conditions and scaffold Young’s modulus values. For each combination of these variables, the explicit equation of the porosity distribution law–i.e the law that describes the pore dimensions in function of the spatial coordinates–was determined that allows the highest amounts of bone to be generated. The results show that the loading conditions affect significantly the optimal porosity distribution. For a pure compression loading, it was found that the pore dimensions are almost constant throughout the entire scaffold and using a FGS allows the formation of amounts of bone slightly larger than those obtainable with a homogeneous porosity scaffold. For a pure shear loading, instead, FGSs allow to significantly increase the bone formation compared to a homogeneous porosity scaffolds. Although experimental data is still necessary to properly relate the mechanical/biological environment to the scaffold microstructure, this model represents an important step towards optimizing geometry of functionally graded scaffolds based on mechanobiological criteria.


Introduction
Functionally Graded Scaffolds (FGSs) for bone tissue engineering are porous biomaterials where the porosity changes with a specific gradient in space. The gradation of porosity enables FGSs to combine together the best mechanical properties of the denser material with those of the more porous one and the resulting material exhibits performances higher than those of the single constitutive materials. Low porosity regions offer high mechanical strength, high porosity regions promote, instead, cell adhesion and support cell growth, proliferation and differentiation [1][2].
Such scaffolds have been successfully utilized in the most variegated domains including the repair of long bone [1,3] and osteochondral [4][5] defects, the maxillofacial [6][7] and the spinal [8] surgery, the cranial reconstruction [9] and the drug delivery systems [1,10]. A large number of studies [11][12][13] are reported in the literature on the manufacturing processes that can be adopted to fabricate these biomaterials. Among the others, the strategy based on the integration of additive manufacturing or rapid prototyping techniques with computer-aided design models seems to be one of the most efficient [2,14]. The possibility of building any scaffold architecture with any type of porosity gradation and the experimental evidence that the geometry of porous scaffolds significantly influences the cellular response and the rate of bone tissue regeneration [15][16][17] led research community to find the possible models that relate the scaffold gradient to the mechanical and biological requirements for the regeneration of the bony tissue [2]. However, to date such models have not been developed yet.
In this article, we attempt to bridge the gap and propose a mechanobiology-driven optimization algorithm that, based on the boundary and loading conditions acting on the scaffold, identifies the best porosity distribution that allows the bone formation to be maximized. Other studies reported in the literature utilized optimization techniques to determine the best scaffold geometry [18][19][20][21][22][23] but none of them adopted mechanobiological criteria and determined the optimal porosity gradient in FGSs. In a previous study [24], the algorithm was utilized to determine the optimal pore dimension in regular structured open-porous scaffolds with homogeneous porosity. In the present study, the model was further developed to include a functionally graded porosity. In particular, three different variables have been investigated: the porosity distribution law, the loading conditions and the scaffold Young's modulus; for each combination of the three variables, the algorithm determines the explicit equation of the porosity distribution law (i.e. the law that describes the pore dimensions in function of the spatial coordinate), that allows the largest volume of the scaffold to be occupied by bone.

Parametric model of an open-porous functionally graded scaffold
The parametric finite element model of an open-porous functionally graded scaffold was created in ABAQUS CAE 1 Version 6.12 (Dassault Systèmes, France). The model has a prismatic geometry with a square t × t = 2548 μm × 2548 μm base and a h = 3822 μm height. The scaffold (represented in yellow, Fig 1A) includes circular pores with a parametric radius A that was assumed to change only along the y direction and remain constant along x and z direction ( Fig  1A and 1B). According to Byrne et al. [25], the scaffold pores were hypothesized to be occupied by granulation tissue (represented in red, Fig 1B). The finite element mesh includes tetrahedral biphasic poro-elastic elements. 4-node linear coupled pore pressure elements (C3D4P) available in ABAQUS were utilized to model both, the scaffold ( Fig 1C) and the granulation tissue ( Fig 1D). The approximate element size was fixed equal to 40 μm.  3,4), pore radius at specific y locations; y max , y min , y int , y int1 , y int2 , specific y locations where the pore radius was determined; m, m i (i = 1, 2, 3), gradient of porosity distribution laws; S, biophysical stimulus regulating the differentiation process; ɣ, octahedral shear strain; ʋ, interstitial fluid flow; ε I , ε II , ε III , principal strains; a, empirical constant a = 3.75%; b, empirical constant b = 3μms -1 ; n resorb , n mature , c, boundaries of the mechano-regulation diagram; A lower , lower bound of the pore radius variability range; A upper , upper bound of the pore radius variability range; V i_bone , volume of the generic element where the formation of mature bone is predicted to take place; n b , number of elements where the formation of mature bone is predicted to take place; V BONE , total volume of the elements where the formation of mature bone is predicted to take place; BO % , percentage of scaffold volume that is occupied by bone; Ω, objective function to optimize; PVPD, Percent Variation of the Pore Dimension; A H , A L , highest and lowest value of A along the y-axis, respectively; BO %_tri-linear , percentage of volume occupied by bone predicted for the tri-linear porosity distribution law; BO %_constant , percentage of volume occupied by bone predicted for the constant porosity distribution law; iBO % , increment of BO % .
Material properties implemented in the finite element model of the granulation tissue are the same as those utilized in previous studies [24,[26][27]. In detail, the Young's modulus was set equal to 0.2 MPa; the permeability to 1×10 −14 m 4 /N/s; the Poisson's ratio to 0.167; the porosity to 0.8; the bulk modulus grain to 2300 MPa; the bulk modulus fluid to 2300 MPa. In order to evaluate the effect of the scaffold mechanical properties on the optimal porosity distribution, three different values of the Young's modulus E were hypothesized: 500, 1000 and 1500 MPa which are the same as those utilized in a previous study [24].
The nodes of the bottom surface of the model were clamped (Fig 1E, 1F and 1G) while those of the upper surface were tied to a rigid plate (represented in blue, Fig 1E, 1F and 1G). For the outer nodes of the granulation tissue the pore pressure was fixed equal to 0 MPa which indicates that the liquid can freely exudate while applying the load. Three different loading conditions were hypothesized in the study: (a) a compression force F V producing a vertical distributed load of F V / (t × t) = 1 MPa ( Fig 1E); (b) a shear force F H producing an horizontal distributed load of F H / (t × t) = 0.5 MPa ( Fig 1F); (c) a mixed compression-shear force F M given by the sum F M ¼ F V þ F H (Fig 1G). The choice of setting F H = 0.5 × F V was done because scaffolds are primarily designed to undergo to compression loading [25]. In all the hypothesized loading conditions, force was ramped over a time period of 1 s that is the possible time in which, a human body motion (such as to assume the erect position or to perform any motion of anatomical regions where a FGS can be implanted), can be completed. The same time interval was utilized in previous studies [24,[28][29].

Porosity distribution laws
The dimension of the circular pores was controlled by the parametric radius A (Fig 1) that was hypothesized to change along the y-direction according to different porosity distribution laws. The coefficients of these distribution laws and hence their gradients were determined via the optimization algorithm described below. The porosity distribution laws considered in the study are the following: constant, linear, bi-linear and tri-linear.
• Constant law. All the pores have the same dimensions (Fig 2A). In this case, the optimization algorithm determines just one coefficient, i.e. A 1 , that is the pore radius of all the scaffold pores.
• Linear law. The dimensions of pore change linearly with y. Two coefficients have to be determined by the optimization algorithm: A 1 and A 2 that are the pore radii at y = y min and y = y max , respectively ( Fig 2B).
• Bi-linear law. The pore radius A changes in the ranges [y min y int ] and [y int y max ] with two different linear laws that assume the same value for y = y int . The coefficients to optimize are three: A 1 , A 2 and A 3 ( Fig 2C).
• Tri-linear law. The dimensions of the pore change in the intervals [y min y int1 ], [y int1 y int2 ], [y int2 y max ] with three different linear laws. The laws defined in the first and second and those defined in the second and third interval assume the same value for y = y int1 and for y = y int2 , respectively. In this case, the optimization algorithm determines four coefficients: A 1 , A 2 , A 3 and A 4 ( Fig 2D).
The specific values of y max , y min , y int , y int1 and y int2 are reported in Fig 2. It is worthy to note that, once the coefficients A i (i = 1, 2, 3, 4) have been determined, the explicit equation of the best porosity distribution, i.e. the equation that describes how the pore radius A changes with y, can be obtained by simply implementing the obtained coefficients in the relationships reported in Table 1.

Computational mechano-regulation model
Once the mesenchymal stem cells invade the scaffold and spread through its pores, the bone regeneration process starts. After dispersal, cells will differentiate. The biophysical stimulus S that regulates the differentiation process was hypothesized to be a function of the octahedral shear strain ɣ and interstitial fluid flow ʋ in the extracellular environment of the cells. In detail, let ε I , ε II , and ε III be the principal strains, the octahedral shear strain ɣ can be defined as: Calling a and b two empirical constants defined as in Huiskes et al. [30], and given by a = 3.75% and b = 3 μms -1 , the biophysical stimulus S can be expressed, according to Prendergast et al. [31], as: Mesenchymal stem cells differentiate into different cell phenotypes according to the following inequalities: where n resorb = 0.01, n mature = 0.53 and c = 3 represent boundaries of the mechano-regulation diagram the values of which are the same as those utilized in other studies [28,[32][33].

Optimization algorithm
The FGS parametric finite element model, the computational mechano-regulation model above described and a numerical optimization routine were combined together in an algorithm written in Matlab 1 (v. R2011b) (Fig 3) that aims to determine, for each of the hypothesized scaffold Young's moduli, loading conditions and porosity distribution laws, the equations of the best porosity distribution that allows the bone formation to be maximized. Considering that 3 scaffold Young's modulus values (i.e. 500, 1000 and 1500 MPa), 3 loading conditions (i.e. F V ,F H , and F M ) and 4 porosity distribution laws (i.e. constant, linear, bi-linear and tri-linear) have been hypothesized, it follows that a total of 3 × 3 × 4 = 36 optimization analyses have been performed in the study.
As a first step, the algorithm requires to select (Block [1]) one of the porosity distribution laws (Block [2]). The initialization of coefficients A i follows (Block [3]), the user can assign to A i initial values that fall within the interval [A lower A upper ], where A lower = 5 μm and A upper = 300 μm have been taken the same as those utilized in a previous study [24]. The algorithm implements the specified initial values of A i into a PYTHON script (Block [4]) that is given in input to ABAQUS. The PYTHON script, based on the values A i , defines in function of the coordinate location y the Table 1. Porosity distribution laws implemented in the study.

Porosity distribution Coefficients to optimize Equation Gradient
Constant law dimension of each pore. The module ABAQUS CAE builds the CAD model of the functionally graded scaffold (Block [5]) with the computed pore dimensions, and after applying the boundary and (one of) the (three) loading conditions above described, generates the finite element mesh (Block [6]). The finite element analysis follows that accounts for geometrical and material nonlinearities (Block [6]). For each element occupying the scaffold pores, i.e. the elements represented in red in Fig 1D, ABAQUS prints (Block [7]) the values of the principal strains ε I , ε II and ε III and of the interstitial fluid flow ʋ that the algorithm utilizes to compute, through the eqs (1) and (2), the magnitude of the biophysical stimulus S (Block [8]). Then, the relationships eq (3) are implemented and for those elements for which the inequality is satisfied, i.e. for those elements where the formation of mature bone is predicted to take place, the volume V i_bone is stored (Block [9]). If n b is the number of elements where inequality eq (4) is satisfied, the algorithm calculates the total volume of these elements as: If V TOT is the total volume of the scaffold model V TOT = t × t × h = 2548 μm × 2548 μm × 3822 μm = 24.814 mm 3 , the algorithm determines the percentage of scaffold volume BO % that is occupied by bone as (Block [10]): and calculates the value of the objective function O as (Block [11]): At this point, the algorithm formulates an optimization problem that includes the coefficients A i as design variables and that aims to minimize the value of the objective function O or, equivalently, to maximize the percentage BO % of volume occupied by bone. It can be claimed, in fact, that the greater the efficiency of the scaffold, the larger the amount of bone produced by the scaffold itself. In an ideal scaffold, 100% of its volume is occupied by bone. The inverse problem described with the eqs (6) and (7) was solved with the Sequential Quadratic Programming (SQP) method available in Matlab, an iterative method for nonlinear optimization. The number of iterations performed by the method can be controlled by means of specific stopping criteria that can be selected by the user and that include a number of tolerances. As one of these stopping criteria is meet, the optimization process ends and after implementing the optimal coefficients A i in the relationships of Table 1, the optimal porosity distributions are traced in function of y (Blocks [14] and [15]). If no stopping criteria are satisfied, the optimization algorithm assigns new values to A i thus generating new candidate solutions (Block [13]). The optimization process terminates when one of the selected stopping criteria is satisfied (Block [12]).
The biophysical stimulus S on which the objective function O depends, was computed based on the hypothesis that the dispersal of mesenchymal stem cells has already taken place and that the only granulation tissue, with the mechanical properties above described, occupies the scaffold pores.
All the computations were performed on a HP Z620-Intel 1 Xeon 1 Processor E5-2620-16Gb RAM. The most expensive optimization analyses were those implementing the tri-linear law that took around 300 hours of computations.

Results
In the case of the compression loading F V the predicted pore dimension experiences small changes (Fig 4A, 4C and 4E) along the y-axis and is almost constant. Independently from the scaffold Young's modulus E, A does not change by more than 15 μm. The general trend (with the exception for the porosity distribution obtained implementing the constant law) that can be observed is that the pore radius in the vicinity of the clamps (i.e. for high values of y) and of the load (i.e. for small values of y) slightly decreases. For increasing values of E, the pore radius, on average, increases. For instance, in the case of E = 500 MPa, the average pore radius is about 190 μm, for E = 1500 MPa, instead, becomes about 220 μm. The percentage of volume occupied by bone BO % increases as we move from the constant to the tri-linear porosity distribution ( Fig  4B, 4D and 4F). Furthermore, increasing values of BO % were predicted for increasing values of the scaffold Young's modulus (Fig 4B, 4D and 4F).
More interesting appears the porosity distribution predicted by the algorithm in the case of the shear load F H (Fig 5) where important changes of the pore dimensions are predicted along the y-axis (Fig 5A, 5C and 5E). The highest values of A are predicted in the vicinity of the load (i.e. for small values of y) while the pore dimensions tend to decrease as we move towards the clamped region. Also in this case BO % increases as we move from the constant to the tri-linear porosity distribution, however, the change of BO % is more significant than in the case of compression load. For increasing levels of E, the average value of BO % increases too (Fig 5B, 5D and 5F).
In the case of mixed load F M , the pore radius A experiences changes that are less important than those predicted in the case of shear load F H but that are certainly larger than those computed in the case of compression load F V (Fig 6A, 6C and 6E). As in the previous case, the pore dimension decreases for increasing values of y. BO % increases as we move from the constant to the tri-linear law and its average value increases for increasing values of the scaffold Young's modulus E (Fig 6B, 6D and 6F). For a fixed value of E and porosity distribution law, the values of BO % predicted in the case of mixed load F M are smaller than those predicted for the other hypothesized loading conditions (Figs 4B, 4D, 4F, 5B, 5D, 5F, 6B, 6D and 6F).
In order to quantify (i) the change of the pore dimensions with y and (ii) the "usefulness" of utilizing a functionally graded scaffold instead of a scaffold with a homogenous porosity distribution we introduced two parameters. The first one, denoted as PVPD, represents the Percent Variation of the Pore Dimension and is defined as: where A H and A L are the highest and the lowest value of A along the y-axis, respectively ( Fig  7A). In general, the higher the PVPD, the larger are the changes of the pore dimension A. The highest values of PVPD have been found in the case of the shear loading F H (Fig 7) where changes of A also by more than 25-30% were predicted (Fig 7C). Slightly lower are the values of PVPD found in the case of the mixed load F M (Fig 7D) and yet less significant those computed in the case of the compression load F V (Fig 7B). Averagely, it appears that PVPD does not depend neither on the scaffold Young's modulus E, nor on the porosity distribution law but does depend on the loading conditions. For the constant law, regardless of the type of load considered, the value of PVPD is zero and is not shown in Fig 7. In general, it appears that as we move from the constant to the linear, bi-linear and, finally, tri-linear porosity distribution law the percentage of volume occupied by bone BO % increases (Figs 4B, 4D and 4F, 5B, 5D and 5F, 6B, 6D and 6F). In particular, the highest values of BO % have been found for the tri-linear law while the lowest ones for the constant law. Therefore, it    Optimization of Functionally Graded Scaffolds makes sense to introduce the second parameter, denoted as iBO % and defined as the increment of BO % when we move from the constant to the tri-linear law. If BO %_tri-linear is the percentage of volume occupied by bone predicted for the tri-linear porosity distribution and BO %_constant the percentage predicted with the constant one, iBO % can be expressed as: As is clear, the higher the values of iBO % , the more "useful" is the utilization of a functionally graded scaffold instead of a homogeneous porosity scaffold. In the limit case where iBO % = 0%, the use of a FGS does not make sense and a homogeneous porosity scaffold has the same potentialities of generating bone as the FG one. On average, the highest values of iBO % were computed in the case of shear loading F H followed by the mixed load F M and the compression load F V , respectively (Fig 8). In particular, among the hypothesized scaffold Young's moduli, the highest values of iBO % were predicted for E = 1000 MPa.
A three-dimensional view of the optimal scaffold geometry predicted for the tri-linear porosity distribution (that is the law with which the highest values of BO % have been obtained) and the shear loading F H is shown in Fig 9. As it can be seen, the pore dimensions change significantly along the y-axis and, on average, increase for increasing values of E.

Discussion
This article presented an optimization algorithm based on mechanobiological criteria and aimed to determine the best porosity distribution in functionally graded scaffolds for bone tissue engineering.
Four porosity distribution laws, three loading conditions and three scaffold Young's moduli were hypothesized. For each combination of these three variables, the optimal microstructure geometry was determined. It was shown that all these variables have a critical effect on the amounts of bone predicted to form within the scaffold pores.
Regarding the porosity distribution law, it was found that designing FGSs with a tri-linear law allows the largest amounts of bone to be generated (Figs 4-6) compared to bi-linear, linear and constant laws. In general, the use of porosity distribution laws with increasing complexity level (i.e. with increasing number of coefficients A i ) leads the scaffold geometry to be better tailored to the specific boundary and loading conditions acting on the construct thus allowing the bone formation to be maximized. Increasing the complexity level of a porosity distribution means, in other words, to include a larger number of design variables and hence, to increase the probability that the optimizer will find a geometry that allows larger amounts of bone to be generated.
More critical appears the effect of the loading conditions. For a pure compression loading, the changes of the pore dimension A are marginal (Figs 4, 7 and 9) and using a FGS allows the formation of amounts of bone slightly larger than those obtainable with a homogeneous porosity scaffold (Fig 8). For a pure shear loading, instead, FGSs allow to significantly increase the bone formation compared to a homogeneous porosity scaffolds (Figs 5B, 5D, 5F and 8) and the pore dimensions change (vs. y) also by more than 20-25% (Figs 7 and 9). This behavior can be justified with the following argument. In the case of pure shear loading, strains increase as we move from the loaded towards the clamped region and hence, the stimulus S, that is a function of the strain, changes in the same manner. In order to maximize the number of elements for which inequality eq (4) is satisfied, the optimization solver tends to reduce the dimension of the pores subjected to higher strain and increase that of the pores subjected to lower strain. In the case of pure compression, instead, (from the macroscopic point of view) the scaffold model is subjected to an uniaxial stress state (with the only exception of the regions close to the loaded Optimization of Functionally Graded Scaffolds and the clamped surfaces where the stress state becomes tri-axial) and then to a more or less uniform distribution of the stimulus S, which explains the approximately uniform dimensions of the pores. The mixed load F M leads to an intermediate situation between the pure compression and the pure shear. Changes of A as well increments of BO % are more important than those predicted in the case of compression force F V but less relevant than those computed with the shear load F H (Figs 6-8).
Finally, regarding the scaffold Young's modulus it appears that the average pore dimension A increases for increasing values of E (Fig 9). This can be justified with the argument that as the Young's modulus increases, the global scaffold stiffness increase too and the optimizer tends to increase the dimensions of the pores to include larger amounts of bone.
To determine the optimal porosity distribution in FGSs some assumptions were made. First of all, the temporal variable was neglected. It was assumed that the scaffold pores are occupied only by granulation tissue, the processes of diffusion of the mesenchymal stem cells and of tissue differentiation were not simulated and the optimization of the porosity distribution was carried out based on the values of the biophysical stimulus registered at the initial time instant. Furthermore, the algorithm does not include scaffold resorption potential [25].
Including the time variable would certainly allow to carry out more accurate predictions on the best porosity distribution but would lead to a dramatic increase of the computational time thus making the algorithm practically not implementable in a "clinical" context. Other aspects such as angiogenesis [34][35][36] and growth factors [37] involved in the process of bone regeneration were not modelled. This model neglects the effect of loads during the initial development of a tissue on a scaffold, i.e. during the phase in which cell attach to the scaffold surface. The scaffold surface is a 2D environment while the model utilized in this study is based on volumetric strains. A model to predict the effect of mechanical signals on cells seeded on the surface of a scaffold has been reported [38]. Another limitation of the model is that a deterministic approach was adopted to determine the biophysical stimulus S, -on the definition of which the optimal porosity distribution law is calculated,-which neglects any possible genetic Optimization of Functionally Graded Scaffolds variability in animal populations. A more general and complete approach would be the probabilistic one and would take into account this variability.
However, despite these limitations, the predictions of the model are consistent with the results of experimental studies. For instance, the patterns of bony tissue predicted in the case of a pure compression load, constant porosity distribution, E = 1000 MPa, are consistent with those of new tissue generated in circular matrix channels observed in histological analyses [39]. In vitro, it was found that, bone forms from the channel walls and tends to growth towards the center of the pore. This same behavior was observed in the numerical model (Fig 10). The grey Fig 10. Patterns of bone predicted in the case of: (i) compression loading; (ii) scaffold Young's modulus E = 1000 MPa; (iii) porosity distribution law: constant. Elements in gray are representative of the regions within the scaffold pores where the algorithm predicts bone formation. Interestingly, the predicted bony tissue patterns appear consistent with those of new tissue formed in three-dimensional matrix channels observed in an in vitro study [39]. Bone formation starts from the pore walls and propagates towards the pore center. Optimization of Functionally Graded Scaffolds elements shown in Fig 10 represent the volumes of the model where the mechano-regulation model predicts the formation of bone. Furthermore, as demonstrated in previous studies a minimum pore size of about 100 μm is required to guarantee a successful bone regeneration process in scaffolds [40]. The pore dimensions predicted by the model are all above this threshold value and well fall within the range of the typical dimensions of the pores of scaffolds for bone tissue engineering [41]. Other studies report that the rate of bone regeneration in scaffold is a function of the scaffold mechanical properties [42]. This is also consistent with the predictions of the present model where the amounts of bone BO % change for changing values of the scaffold Young's modulus (Figs 4-6 and 9).

Conclusions
A mechanobiology-driven optimization algorithm was presented to determine the optimal porosity distribution in functionally graded scaffolds. The results presented in this paper show that the loading conditions are pivotal in determining optimal porosity distribution. For a pure compression loading, it was predicted that the changes of the pore dimension are marginal and using a FGS allows the formation of amounts of bone slightly larger than those obtainable with a homogeneous porosity scaffold. For a pure shear loading, instead, FGSs allow to significantly increase the bone formation compared to a homogeneous porosity scaffold. Increasing pore dimensions are predicted for increasing values of the scaffold Young's modulus. Increasing the number of coefficients that define a porosity distribution law allows to design more performing scaffolds capable of generating larger amounts of bone.
The model predictions appear reasonably consistent with what is observed in vitro. Although experimental data is still necessary to properly relate the mechanical/biological environment to the scaffold microstructure geometry, this model represents an important step towards optimizing geometry of functionally graded scaffolds and/or stimulation regimes based on mechanobiological criteria.