A mechanobiological computer optimization framework to design scaffolds to enhance bone regeneration

The treatment of large bone defects is a clinical challenge. 3D printed scaffolds are a promising treatment option for such critical-size defects. However, the design of scaffolds to treat such defects is challenging due to the large number of variables impacting bone regeneration; material stiffness, architecture or equivalent scaffold stiffness—due it specific architecture—have all been demonstrated to impact cell behavior and regeneration outcome. Computer design optimization is a powerful tool to find optimal design solutions within a large parameter space for given anatomical constraints. Following this approach, scaffold structures have been optimized to avoid mechanical failure while providing beneficial mechanical stimulation for bone formation within the scaffold pores immediately after implantation. However, due to the dynamics of the bone regeneration process, the mechanical conditions do change from immediately after surgery throughout healing, thus influencing the regeneration process. Therefore, we propose a computer framework to optimize scaffold designs that allows to promote the final bone regeneration outcome. The framework combines a previously developed and validated mechanobiological bone regeneration computer model, a surrogate model for bone healing outcome and an optimization algorithm to optimize scaffold design based on the level of regenerated bone volume. The capability of the framework is verified by optimization of a cylindrical scaffold for the treatment of a critical-size tibia defect, using a clinically relevant large animal model. The combined framework allowed to predict the long-term healing outcome. Such novel approach opens up new opportunities for sustainable strategies in scaffold designs of bone regeneration.

1 Introduction

Large bone defects do not heal spontaneously if they exceed a critical dimension. The current gold standard treatment in such clinical cases consists of autologous bone grafting but has several drawbacks such as the need for a second surgery and a limited amount of adequate bone graft tissue (Dimitriou et al., 2011; Schlundt et al., 2018). 3D-printed scaffolds appear as a promising alternative treatment strategy for large bone defects, with good pre-clinical results (Reichert et al., 2012; Pobloth et al., 2018). However, their translation to the clinic has been limited so far, partly due to the large number of scaffold design variables influencing the regeneration process, which limits the predictability of the clinical outcome; hence, bone scaffold design is often the result of a trial-and-error process.

Several research groups have departed from a trial-and-error approach by using computational optimization methods when designing bone scaffolds. A common approach consisted in maximizing or minimizing specific mechanical properties such as stiffness to mimic the material properties of bone tissue (Almeida and da Silva Bártolo, 2010; Chen et al., 2011; Xiao et al., 2012; Dias et al., 2014; Wang et al., 2016; Metz et al., 2020). Such concept may be considered limited, since endogenous bone regeneration is limited in too stiff environments (Prendergast et al., 1997; Pobloth et al., 2018). Other groups have targeted specific values for stiffness and/or diffusivity of the scaffold to mimic the tissue being replaced (Hollister et al., 2002; Hollister and Lin, 2007; Sturm et al., 2010; Wieding et al., 2014; Makowski and Kuś, 2016; Chang et al., 2017). Nonetheless, scaffold properties mimicking a missing tissue do not ensure being supportive to the endogenous bone regeneration process (Petersen et al., 2018).

Few studies have focused on the mechanobiological potential of a scaffold design as a means of optimization. Boccaccio and colleagues optimized periodic scaffolds to provide favorable mechanical conditions for bone regeneration immediately after implantation (Boccaccio et al., 2016a, 2016b, 2018a, 2018b; Percoco et al., 2020). However, such approach is also limited as the bone regeneration process is known to be highly dynamic, where initial optimal conditions cannot ensure an optimal regeneration outcome (Perier-Metz et al., 2021). Over the course of bone healing, different types of tissues are formed, remodeled and resorbed in specific areas of a defect, creating a highly dynamic local mechanical niche that changes over time and further influences cell response and tissue formation. A few studies confirmed that designs optimized for the situation immediately after surgery did not yield optimal bone regeneration in various set-ups (Bashkuev et al., 2015; Perier-Metz et al., 2021; Wu et al., 2021). These findings suggest that bone scaffold optimization should rather be performed on an intended regeneration outcome and not only the initial post-surgery setting, preferably as a time-dependent mechanobiological optimization (Metz et al., 2020).

So far, to the best of our knowledge, only two time-dependent, mechanobiology-based bone scaffold optimization studies have been performed. In the first one, Poh et al. (2019) optimized a scaffold porosity distribution for a large defect in a sheep long bone, based on a simplified bone regeneration model. However, their model did not include any other tissue type than bone (e.g., cartilage, fibrous tissue). Furthermore, the optimization set-up was restricted to 1D, meaning that the actual scaffold micro-structure was not optimized. More recently, Wu and colleagues proposed a method for a time-dependent mechanobiology-based topology optimization of bone scaffolds (Wu et al., 2021). However, also this model did not account for any other tissue type but bone nor for the needed initial cell infiltration into the defect.

To overcome these limitations, we propose a novel in silico framework to optimize scaffold design based on the level of regenerated bone volume aimed at the end of the regeneration process. Such prediction has to consider the activity of individual cells and the potential for different tissue types to be formed during the healing process. To achieve this, a computer framework combining a multiscale mechanobiological bone regeneration model, a surrogate model for bone healing outcome and an optimization algorithm was developed. This framework was applied to optimize scaffolds for a large defect in a sheep tibia.

2 Materials and methods

To realize an optimization taking into account the regeneration outcome, this study focuses on the design optimization of a cylindrical scaffold to be implanted in a large bone defect in the sheep tibia as an example. The in silico framework to optimize scaffolds was realized by using a finite element model that was coupled to an agent-based model to simulate bone regeneration within such scaffolds. This computer model was first used to simulate bone regeneration in a relatively large set of scaffold designs, to generate the interpolation data necessary for a surrogate model. Based on the results of the bone regeneration model, a simplified relationship between the scaffold design parameters and the predicted regenerated bone outcome was computed (i.e., surrogate model). Thereafter, this relationship was used in an optimization framework, to find the optimal scaffold design parameters to maximize the amount of regenerated bone.

2.1 Finite element model

A finite element model was developed where the geometry of the defect and the surrounding bone extremities was based on previously published in silico-in vivo studies (Pobloth et al., 2018; Perier-Metz et al., 2020). A 4-cm segmental defect in a sheep tibia was modelled by a cylinder of outer radius 10 mm (cortical bone) and inner radius 7.5 mm, containing the bone marrow (Figure 1A). The defect was modelled to be stabilized by an internal stainless steel fixation plate and six screws, following previous experimental setups (Pobloth et al., 2018). A regenerating zone, the callus, was modelled in and around the defect by rotating a circle arc of maximum width 10 mm at mid-height and overlapping the intact bone extremities by 10 mm.

FIGURE 1

FIGURE 1. (A) Finite element model set-up: longitudinal cut through the defect configuration showing the intact bone extremities, the implanted scaffold and the fixation plate; the dashed line represents the symmetry plane used in the analysis. (B–D) Cylindrical scaffold description and parametrization: (B) longitudinal view of the full scaffold geometry; (C) radial view of the scaffold defining the vertical pore size parameter x₃; (D) longitudinal view of the half scaffold (taking advantage of the symmetry) defining the horizontal pore size parameters x₁ and x₂.

A scaffold with outer radius 10 mm and inner radius 5 mm and pores with square sections following three perpendicular directions was modelled (Figures 1B,C). The pore centers had a spacing of 1 mm in both horizontal directions (longitudinal view, Figures 1B,D); whereas the 12 pores extruded along the vertical direction were placed regularly in a circle (radial view, Figure 1C). Three variables were defined to parametrize the scaffold geometry: x₁ was the pore size closest to the bone extremities and x₂ the pore size at mid-height of the scaffold, assuming a linear gradient of the pore sizes in the three rows between these (Figure 1D); and x₃ defined the pore size in the radial plane (Figure 1C). The parameters x₁ and x₂ were constrained between 0.3 and 3.8 mm and x₃ between 0.3 and 3 mm. The lower bound was dictated by biological requirements to ensure good vascularization and nutrient supply, as experimental studies suggested a minimal pore size of 300 µm (De Witte et al., 2018), while the upper bound was chosen to avoid that neighboring pores overlap.

This model was implemented in the finite element (FE) software Abaqus v.6.18 (Simulia, Rhode Island) to perform mechanical analyses. Because of the symmetry of the geometry, only one half of the model was simulated, following the dashed line along the horizontal plane at mid-height (Figure 1A). The different geometrical parts were meshed with second-order elements as follows: hexahedral elements of average size 2.5 mm for the intact cortical bone, the bone marrow and the plate; beam elements of size 1 mm for the screws; and tetrahedral elements of average sizes 1.2 mm for the callus and scaffold—this size was adapted depending on the size of the meshed geometry, in particular the scaffold struts.

The following constraints were defined between the different geometrical parts: the intact bone extremities (cortical bone and bone marrow) were tied to the callus; the screws were tied to the intact bone; the screw heads were constrained by multi-point constraints of type beam with the plate holes; the callus was constrained such that it was tied to the scaffold.

The biological parts (cortical bone and bone marrow) were modelled as poroelastic materials with properties summarized in Table 1. The callus space was assumed to be initially filled with granulation tissue (poroelastic material, Table 1). Over time, regenerating tissue material properties were updated according to the predicted tissue formation of the bone regeneration model (Section 2.2).

TABLE 1

TABLE 1. Material properties (Perier-Metz et al., 2020).

A compression load of 1372 N and a bending moment of 17.125 Nm were applied on the proximal end of the bone, corresponding to two body weights (BW) proximal-distal compression and an anterior-posterior moment of 0.025 BWm at the fixed end of a 20-cm intact bone, respectively (Duda et al., 1997; Perier-Metz et al., 2020). A symmetry boundary condition was defined on the xy symmetry plane (“ZSYMM”). Lastly, pore pressure was constrained to be zero on the outer surfaces of the poroelastic parts (callus, cortical bone and bone marrow).

Two different scaffold optimization studies were performed, assuming different scaffold Young’s moduli (Table 1):

1) titanium, defined as a linear elastic material with Young’s modulus 104 GPa

2) a hypothetical very soft material, defined as a linear elastic material with Young’s modulus 0.2 MPa (similar to the Young’s modulus of granulation tissue)

2.2 Bone regeneration model

A previously described mechanobiological bone regeneration (MBBR) model, which was able to explain bone regeneration within scaffolds in different experimental setups, was used to predict tissue formation during the healing process (Perier-Metz et al., 2020, 2022). It consisted in a 3D agent-based computer model implemented in C++ and coupled with the afore-mentioned FE model (Section 2.1). Each agent, of size 100 μm, could represent one cell of one of the following phenotypes: progenitor cell, fibroblast, chondrocyte and immature or mature osteoblast, with the associated extracellular matrix: granulation tissue, fibrous tissue, cartilage and immature or mature bone, respectively.

Progenitor cells were allowed to migrate randomly with an average speed of 30 μm/h (Appeddu and Shur, 1994). The origin of these cells was mimicked by initially seeding them randomly in the bone marrow cavity and along the periosteum with a 30% occupancy rate (Perier-Metz et al., 2020). The rest of the tissue volume was considered cell-free at the initial time point.

Differentiation, proliferation and apoptosis of the cells depended on a mechanoregulation stimulus S based on octahedral shear strain $\gamma$ and fluid velocity $\nu$, defined by $S=\frac{\gamma }{a}+\frac{\nu }{b}$ , with the empirical parameters a and b: $a=0.0375$ and $b=0.03mm.s^{-1}$ (Huiskes et al., 1997; Lacroix and Prendergast, 2002). Mechanical thresholds were defined for the stimulus (S) favoring bone resorption (S < 0.01), mature bone (0.01 < S < 0.53), immature bone (0.53 < S < 1), cartilage (1 < S < 3) or fibrous tissue (S > 3) (Perier-Metz et al., 2020). Cells would proliferate with a predefined proliferation rate (60%, 55%, 20% or 30% for progenitor cells, fibroblasts, chondrocytes and osteoblasts, respectively) providing that they were in the right mechanoregulation stimulus; or would undergo apoptosis otherwise, with a different rate (5%, 5%, 10% or 16% for progenitor cells, fibroblasts, chondrocytes and osteoblasts, respectively) (Perier-Metz et al., 2020). Progenitor cells differentiated with a rate of 0.3 days⁻¹ towards the cell phenotype corresponding to the local stimulus. In addition, surface-guided tissue deposition was implemented, meaning that differentiation could happen only following existing surfaces (tissues or scaffold), and graft was assumed to be present in the scaffold pores where it stimulated proliferation and differentiation, as previously described (Perier-Metz et al., 2020).

This agent-based model was coupled with the FE model in two ways: 1) the mechanoregulation stimulus that dictated the cell behaviors was derived element-wise from the FE analysis; and 2) the tissue material properties were updated in the callus of the FE model at every iteration depending on the agent-based model cell distribution. To do so, each element was mapped to the agents it contained and its material properties were defined according to a rule of mixtures: the average material properties of the tissues predicted in the agents within each single element was computed, and the obtained value was further averaged over the last ten iterations to account for the time needed for tissue deposition and maturation (Lacroix and Prendergast, 2002). The FE analysis and agent-based simulations were run iteratively to predict the full regeneration process, where one iteration represented one day.

2.3 Surrogate optimization set-up

The bone regeneration model described above (Section 2.2) is computationally intensive (it takes around 6 h on a standard workstation to simulate one full regeneration process). However, for the optimization process, the regeneration outcome needs to be predicted multiple times, until the optimum scaffold is found. Therefore, to ensure that the optimization process runs in a reasonable time, a surrogate optimization approach was adopted. A surrogate computer model (a simplified input-output relationship) of the bone regeneration model described above (Section 2.2) was built based on a set of initially sampled input parameters obtained by a design of experiments technique. The input parameters were the scaffold design variables (x₁-x₃) whereas the output was defined as the predicted fraction of regenerated bone volume within the scaffold pores after 24 weeks. Thereafter, the scaffold design optimization was performed using the surrogate model. To ensure surrogate model accuracy and optimality of the outcome, the optimum design predicted by the surrogate model was computationally tested for bone regeneration outcome using both the surrogate model and the MBBR model. A first loop of the framework was used to ensure that the difference in the regenerated bone volume value between the surrogate and the MBBR model was below 5% (Figure 2A); as long as the error was bigger, the input values and the corresponding output (as simulated with the MBBR model) were added to the dataset to compute a new surrogate model, and a new optimum was determined based on this improved surrogate model. The second loop consisted in building the surrogate model with the given dataset and computing an optimum based on this surrogate model until an optimum was reached that was better than all previously stored values (to avoid local optima). The output was simulated again with the MBBR model to ensure a good confidence in the result. If a new loop was necessary, this new data point was added to the dataset to build the surrogate model and improve it further.

FIGURE 2

FIGURE 2. (A) Flow chart of the scaffold design optimization computational framework for enhanced bone regeneration. (B) Locations of the initial samples (in the space defined by the three parameters x₁, x₂, x₃) for which the MBBR model is run to compute the surrogate model.

The initial dataset used to build the surrogate model was determined using the Latin Hypercube Sampling technique (Simon et al., 2002). Based on preliminary studies, 20 times the number of variables was found to be sufficient to obtain an accurate enough surrogate model, i.e, 60 samples in this set-up. The locations of the samples are given in Figure 2B.

The surrogate model was built using the predicted regenerated bone volume using the Matlab kriging toolbox DACE (Lophaven et al., 2002). An exponential auto-correlation was assumed as it gave the best results in a preliminary study. To perform the optimization, the Matlab Global optimization toolbox was used where different optimization algorithms were tested: direct search, genetic algorithm and particle swarm optimization (The MathWorks, Inc. 2020). For the direct search, an initial guess for x1, x2 and x3 has to be defined to initialize the optimization process: the values 2.05, 2.05, and 1.65 mm were arbitrarily used. In all algorithms, the default settings were used to perform the optimization.

In summary, the surrogate optimization framework includes the following steps (Figure 2A):

1) Run the MBBR model for a set of input parameters obtained by Latin Hypercube Sampling technique (LHS).

2) Build a kriging surrogate model based on the initial data using the Matlab kriging toolbox DACE.

3) Perform optimization based on the surrogate model predictions using the Matlab Global optimization toolbox (The MathWorks, Inc. 2020): direct search, genetic algorithm (GA), particle swarm optimization (PSO)

4) Run the MBBR model for the optimum scaffold parameters found in step 3 and compare the MBBR model output to the surrogate prediction (loop 1) or current optimal value (loop 2)

5) Repeat steps two to four to ensure model accuracy (loop 1, allowing max 5% discrepancy between the surrogate prediction and the MBBR model output) and optimality of the outcome (loop 2)

To automate the creation of the FE models for each scaffold design, a Python script was developed for Abaqus, which allowed to define the pore corner positions based on the pore size. MATLAB R2018b (The MathWorks, Inc. 2020) was used to run the actual optimization framework, including launching Abaqus CAE for the geometry update and the C++-Abaqus MBBR model. Pore size values (input) and corresponding output (bone volume in the scaffold pores) were stored in a text file for further analysis.

2.4 Analysis of the outcome

The main outcome of the MBBR simulations was the predicted regenerated bone volume fraction in the scaffold pores after 24 weeks. Regenerated bone volume fraction was computed as the total amount of agents within the scaffold pores occupied by osteoblasts divided by the total number of agent positions within the scaffold pores.

The porosity of the resulting scaffold was computed as the scaffold pore volume divided by the hollow cylinder volume (outer radius 10 mm and inner radius 5 mm). The regenerated bone volume fraction could then be related to the scaffold macroscopic porosity.

In addition, histology-like images of bone regeneration predictions were computed in the mid-sagittal plane of the defect using histology-like colors, similar to the experimental Safranin Orange/von Kossa staining: bone in black, fibrous tissue in light red and cartilage in dark red (Pobloth et al., 2018; Perier-Metz et al., 2020). This representation allowed a comparison between different designs and study cases.

3 Results

3.1 Titanium scaffold optimization

The first step of the optimization process of the titanium scaffold consisted in simulating the 60 initial designs defined by the Latin hypercube sampling strategy (Figure 2B).These simulations yielded a regenerated bone volume fraction in the scaffold pores varying from 10% to 95%. The corresponding designs had porosities ranging from 13% to 95%, with a tendency towards a better regeneration potential for higher porosities until ca. 90% (Figure 3A). Above this porosity, the designs yielded very bad healing outcomes due to high deformations stimulating fibrocartilage formation instead of bone. Nevertheless, the porosity alone was not a good indicator of the healing potential, as the same porosity yielded many different healing outcomes. For instance, 50%-porosity scaffolds yielded 35%–90% regenerated bone volume in the scaffold pores (Figure 3A).

FIGURE 3

FIGURE 3. Predicted regenerated bone volume fraction (bone volume/total pore volume) within the scaffold pores for scaffold designs with different porosities for 60 initial scaffold designs: (A) made of titanium; (B) made of a soft material (Young’s modulus 0.2 MPa).

The three tested optimization algorithms (direct search, GA and PSO) suggested similar optimal titanium scaffold designs yielding ca. 96% regenerated bone (Table 2). The best design was obtained with the PSO algorithm, where the scaffold design was defined by $x_1=3.12mm$ , $x_2=3.26mm$ and $x_3=2.6mm$ This scaffold had a porosity of 85% and yielded roughly 96% of regenerated bone in the scaffold pores (Figure 4A). The optimal scaffold did not show a clear advantage of a gradient in the scaffold pore size for bone regeneration, since the values of x₁ and x₂ were very close to each other. Interestingly, one of the very graded designs ($x_1=1.84mm$ , $x_2=3.45mm$ , $x_3=2.26mm$), that was part of the interpolation dataset, performed nearly as well, with 93% regenerated bone but only 71% porosity (Figure 4B). Compared to the optimal design, some areas of bone resorption appeared around some scaffold walls.

TABLE 2

TABLE 2. Optimization results for both scaffold compositions using various algorithms.

FIGURE 4

FIGURE 4. Scaffold design, initial absolute principal strain distribution in the mid-sagittal plane and 24-week histology predictions in the mid-sagittal plane for different scaffold designs: (A) optimal titanium scaffold design; (B) a good titanium graded design (${\mathrm{X}}_1=1.84\mathrm{mm}$ ${\mathrm{X}}_2=3.45\mathrm{mm}$, ${\mathrm{X}}_3=2.26\mathrm{mm}$); (C) scaffold design optimized for titanium but now made of soft material; (D) optimal soft material scaffold design.

3.2 Soft scaffold optimization

The predicted regenerated bone volume fractions for the interpolation dataset using soft material properties for the scaffold varied from 20% to 60%, thus showing a reduced healing potential compared to the titanium scaffolds. The corresponding designs of the interpolation dataset had the same geometries as the ones used for the titanium scaffold optimization; therefore, their porosities ranged from 13% to 95% as well. A tendency towards a lower regeneration potential for higher porosities was observed (Figure 3B). Also in this case, the porosity alone was not a good indicator of the healing potential, as e.g., 50%-porosity scaffolds yielded 15%–50% regenerated bone volume in the scaffold pores (Figure 3B).

The optimal titanium scaffold design was predicted to perform very badly when assuming that it would be made of the very soft material with Young’s modulus 0.2 MPa: large amounts of fibrous tissue and cartilage were predicted to grow. As a consequence, only 27% of bone was present in the scaffold pores after 24 weeks compared to 96% when the scaffold was made of titanium (Figure 4C). The much softer material yielded high strain values in the defect, thus leading to a mechanoregulation stimulus favoring more fibrocartilage formation.

The optimal soft scaffold was obtained with the PSO algorithm and the scaffold design was defined by $x_1=1.9mm$, $x_2=0.77mm$ and $x_3=0.7mm$ (Table 2). This design was more graded than the titanium optimum (larger pores close to the intact bone extremities), but far less porous (24% porosity). It yielded only 67% regenerated bone in the scaffold pores after 24 weeks, but large amounts of fibrocartilage in the center of the defect (Figure 4D). The other two algorithms (direct search, GA) converged to different optima which were found to yield even less regenerated bone (51% and 65%, respectively) at similar levels of porosity (Table 2).

3.3 Computational efficiency

The optimization processes presented in Sections 3.1 and 3.2 took ca. 2 weeks to perform, including 10 days for the simulation of the bone regeneration process in the 60 initial scaffold geometries used to build the surrogate model and 2–7 days for each tested optimization algorithm (5–20 additional MBBR simulations) on a standard desktop PC. In particular, when using the PSO algorithm—the algorithm that performed best in both cases—for the optimization step, 100 h (ca. 4 days) of computation were needed for the titanium scaffold optimization and 26 h for the soft material scaffold, in addition to the initial 10 days for the initial dataset simulations. In comparison, a standard optimization approach would have required thousands of simulation runs lasting 6 h each, i.e., years of computations. Hence, the computing performance was improved by a factor of at least 100.

4 Discussion

We propose here a computational framework for time-dependent mechanobiological optimization of 3D-printed scaffolds towards enhanced bone regeneration. Our method is based on the bone regeneration outcome to optimize scaffold design, instead of only the post-surgery mechanical stimulus or scaffold mechanical properties. The framework uses surrogate modelling and design of experiments methods that allow for a set-up that can be run in a reasonable amount of time on standard computing machines. In this study, the framework was used to optimize a scaffold for the treatment of a large bone defect in sheep, where scaffolds of two different material properties were considered: titanium and a very soft material (with a Young’s modulus similar to that of granulation tissue).

The optimization of the titanium scaffold design suggested a very porous design with large pores (86% porosity) to be optimal for bone regeneration, yielding 96% regenerated bone in the scaffold pores. This result is in good agreement with an in silico study performed by Byrne and colleagues that found higher porosities to be more beneficial for bone regeneration (Byrne et al., 2007). A previous pre-clinical study has also shown that a honeycomb titanium scaffold with 84% porosity leads to bone regeneration within the here simulated large bone defect (Pobloth et al., 2018). For the experimental scaffold, the computer model of bone regeneration previously predicted a 98% regenerated bone in the scaffold pores (Perier-Metz et al., 2020). Even though the optimization was performed on a relatively simple scaffold geometry with square-section pores, the computed optimum proved capable of performing as good as a more complex experimental design. This first result is therefore very promising as the optimum of a simple parameterized geometry is already predicted to yield very good bone regeneration capabilities. Further studies should test the potential of the optimization approach in other, more complex scaffold designs, which could have a larger number of design parameters and therefore yield even better results.

A sensitivity analysis of the optimal titanium scaffold design to the loading conditions was further performed to assess the validity of the optimization process (Supplementary Material). Indeed, the optimal design was shown to still yield mostly bone regeneration under 50% higher or lower loading conditions, what suggested a good robustness of the design under loading uncertainties. Other bone scaffold optimization studies have shown that the optimal scaffold design (according to the initial mechanical stimulus) depended on the assumed loading conditions (Boccaccio et al., 2016a; 2016b; 2018a). Here, the lower loading conditions decreased more strongly the regenerated bone volume fraction within the scaffold pores (79% compared to 96%).

The optimization of the soft scaffold (Young’s modulus of 0.2 MPa) showed a very different result, with the optimum being far less porous (24%) and yielding less bone than the optimal titanium scaffold (67%). Interestingly, the design was graded, with larger pores close to the intact bone extremities and smaller ones in the center of the defect. This result is in line with a previous in silico study optimizing the porosity distribution of a large bone defect scaffold (Poh et al., 2019). A reason for this distribution of pore sizes might be the need to sustain higher strains in the center of the defect, as bone regenerates from the intact extremities towards the center so that the remaining granulation tissue in the center of the scaffold is compressed together from the bone formed above and below. Here, the chosen material was very soft; the experimental literature suggests intermediate material properties as an alternative to foster bone regeneration, e.g., polymers and/or ceramics (Reichert et al., 2012; Reznikov et al., 2019). Using our set-up, future studies could perform a polymer-ceramic composite scaffold optimization for large bone defect regeneration, including the degradable behavior of these materials.

Two previous studies have used mechanobiological algorithms to perform a time-dependent bone scaffold optimization: topology optimization studies in 2D (large defect) and 3D (partial defect) (Wu et al., 2021) and a porosity distribution optimization study in 1D (Poh et al., 2019). In the first case, the use of topology optimization techniques allowed a free design, where no topology was imposed a priori. However, Wu and colleagues did not take into account the actual cell invasion of the scaffold pores nor the different tissue types taking part in the bone regeneration process (e.g., fibrous tissue, cartilage). In addition, they assumed a periodic design and optimized only one (cubic or square) representative volume element of this scaffold, thus not allowing any variations between different scaffold regions (Wu et al., 2021). The second study consisted in a very fast method as it focused only on one dimension; however, it was strongly limited as the optimal scaffold micro-structure could therefore not be predicted (Poh et al., 2019). In fact, this method should be combined with a micro-structure optimization step for a practical usage by means of a multi-scale optimization (Dondl et al., 2019).

This study had some limitations related to the mechanobiological bone regeneration model. Vascularization was not included, although it is known to be a very limiting factor in large bone defect regeneration; however, this limitation was mitigated by the use of large pores to avoid impairing the re-vascularization process (minimum 300 µm) and the assumption that graft was implanted into the scaffold, which is likely to enhance the vascularization process. Moreover, the bone regeneration computer model used here has been shown to have good predictive capabilities in two independent experimental setups where different scaffolds designs were used (Perier-Metz et al., 2020, 2022). Another limitation was that each agent point was bigger than actual cell sizes (Isaksson et al., 2008), so that an agent was assumed to contain several cells and the corresponding tissue; tissue deposition was therefore not implemented as such. Preliminary studies showed that this simplification did not affect much the final regeneration predictions but rather their dynamics, what was not the focus of the present study. The relatively coarse mesh size is a further limitation of the MBBR model, which might in particular result in distorted predictions of the tissue composition and material properties due to the averaging of material properties over a full element; however, preliminary mesh convergence studies performed in a previous study showed that the mesh size did not impact the mechanoregulation stimulus prediction (Perier-Metz et al., 2020). In this study, optimization analyses were done for two specific scaffold materials, with the very stiff titanium on the one hand [comparable to (Pobloth et al., 2018)] and the very soft granulation tissue-like material on the other (Petersen et al., 2018); future studies should take into account degradable polymeric or composite (polymer-ceramic) materials that are considered for clinical usage. This would add a further dynamic aspect in the optimization process due to the degradation of the scaffold material. Also, it would have been possible to include the material stiffness in the optimization process, but this might be less realistic (only some stiffness values are possible) and would distort the optimization results, as the variation of material stiffness values can be of orders of magnitude compared to the pore size values. Lastly, despite the implementation of a surrogate modelling approach, the computation time (a few weeks for one scaffold design optimization set-up) is still too long for a clinical routine usage, where an optimal design should be designed in a few days for a patient-specific case. Future studies should take advantage of high-performance computing technologies by using dedicated hardware and parallelized algorithms to further reduce the computing time.

In summary, we propose here a technological platform that allows to optimize hierarchically-structured bone scaffold designs not only against mechanical failure and initial in-growth of bone but for a sustainable long-term optimized regenerative process. With improved computational efficiency and providing that the healing potential of an individual is known, this method could be employed for the development of personalized 3D-printed bone scaffolds to ensure an optimal regeneration outcome to a given patient.

Data availability statement

The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

Author contributions

All authors contributed to the conception and design of the study. CP-M and SC developed the in silico models. CP-M collected and interpreted the data. CP-M wrote the first draft of the manuscript. SC provided technical guidance on the project. All authors revised and edited the manuscript and approved its content.

Funding

This work was part of Camille Perier-Metz’ PhD project funded by MINES Paris—PSL Research University. The authors acknowledge funding from the German Federal Ministry of Education and Research (BMBF) (grant number: 01ZX 1910).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fbioe.2022.980727/full#supplementary-material

References

Almeida, H. de A., and da Silva Bártolo, P. J. (2010). Virtual topological optimisation of scaffolds for rapid prototyping. Med. Eng. Phys. 32, 775–782. doi:10.1016/j.medengphy.2010.05.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Appeddu, P. A., and Shur, B. D. (1994). Molecular analysis of cell surface beta-1, 4-galactosyltransferase function during cell migration. Proc. Natl. Acad. Sci. U. S. A. 91, 2095–2099. doi:10.1073/pnas.91.6.2095

PubMed Abstract | CrossRef Full Text | Google Scholar

Bashkuev, M., Checa, S., Postigo, S., Duda, G., and Schmidt, H. (2015). Computational analyses of different intervertebral cages for lumbar spinal fusion. J. Biomech. 48, 3274–3282. doi:10.1016/j.jbiomech.2015.06.024

PubMed Abstract | CrossRef Full Text | Google Scholar

Boccaccio, A., Fiorentino, M., Uva, A. E., Laghetti, L. N., and Monno, G. (2018a). Rhombicuboctahedron unit cell based scaffolds for bone regeneration: Geometry optimization with a mechanobiology - driven algorithm. Mater. Sci. Eng. C 83, 51–66. doi:10.1016/j.msec.2017.09.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Boccaccio, A., Uva, A. E., Fiorentino, M., Bevilacqua, V., Pappalettere, C., and Monno, G. (2018b). “A computational approach to the design of scaffolds for bone tissue engineering,” in Advances in bionanomaterials, bionam 2016. Editors S. Piotto, F. Rossi, S. Concilio, E. Reverchon, and G. Cattaneo (Berlin: Springer-Verlag Berlin), 111–117.

CrossRef Full Text | Google Scholar

Boccaccio, A., Uva, A. E., Fiorentino, M., Lamberti, L., and Monno, G. (2016a). A mechanobiology-based algorithm to optimize the microstructure geometry of bone tissue scaffolds. Int. J. Biol. Sci. 12, 1–17. doi:10.7150/ijbs.13158

CrossRef Full Text | Google Scholar

Boccaccio, A., Uva, A. E., Fiorentino, M., Mori, G., and Monno, G. (2016b). Geometry design optimization of functionally graded scaffolds for bone tissue engineering: A mechanobiological approach. PLOS ONE 11, e0146935. doi:10.1371/journal.pone.0146935

PubMed Abstract | CrossRef Full Text | Google Scholar

Byrne, D. P., Lacroix, D., Planell, J. A., Kelly, D. J., and Prendergast, P. J. (2007). Simulation of tissue differentiation in a scaffold as a function of porosity, Young’s modulus and dissolution rate: Application of mechanobiological models in tissue engineering. Biomaterials 28, 5544–5554. doi:10.1016/j.biomaterials.2007.09.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Chang, C.-C., Chen, Y., Zhou, S., Mai, Y.-W., and Li, Q. (2017). “Computational design for scaffold tissue engineering,” in Biomaterials for implants and scaffolds. Editors Q. Li, and Y. W. Mai (Berlin: Springer-Verlag Berlin), 349–369.

CrossRef Full Text | Google Scholar

Chen, Y., Zhou, S., and Li, Q. (2011). Microstructure design of biodegradable scaffold and its effect on tissue regeneration. Biomaterials 32, 5003–5014. doi:10.1016/j.biomaterials.2011.03.064

PubMed Abstract | CrossRef Full Text | Google Scholar

De Witte, T.-M., Fratila-Apachitei, L. E., Zadpoor, A. A., and Peppas, N. A. (2018). Bone tissue engineering via growth factor delivery: From scaffolds to complex matrices. Regen. Biomater. 5, 197–211. doi:10.1093/rb/rby013

PubMed Abstract | CrossRef Full Text | Google Scholar

Dias, M. R., Guedes, J. M., Flanagan, C. L., Hollister, S. J., and Fernandes, P. R. (2014). Optimization of scaffold design for bone tissue engineering: A computational and experimental study. Med. Eng. Phys. 36, 448–457. doi:10.1016/j.medengphy.2014.02.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Dimitriou, R., Mataliotakis, G. I., Angoules, A. G., Kanakaris, N. K., and Giannoudis, P. V. (2011). Complications following autologous bone graft harvesting from the iliac crest and using the ria: A systematic review. Injury 42, S3–S15. doi:10.1016/j.injury.2011.06.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Dondl, P., Poh, P. S. P., Rumpf, M., and Simon, S. (2019). Simultaneous elastic shape optimization for a domain splitting in bone tissue engineering. Proc. R. Soc. A 475, 20180718. doi:10.1098/rspa.2018.0718

CrossRef Full Text | Google Scholar

Duda, G. N., Eckert-Hübner, K., Sokiranski, R., Kreutner, A., Miller, R., and Claes, L. (1997). Analysis of inter-fragmentary movement as a function of musculoskeletal loading conditions in sheep. J. Biomechanics 31, 201–210. doi:10.1016/S0021-9290(97)00127-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Hollister, S. J., and Lin, C. Y. (2007). Computational design of tissue engineering scaffolds. Comput. Methods Appl. Mech. Eng. 196, 2991–2998. doi:10.1016/j.cma.2006.09.023

CrossRef Full Text | Google Scholar

Hollister, S. J., Maddox, R. D., and Taboas, J. M. (2002). Optimal design and fabrication of scaffolds to mimic tissue properties and satisfy biological constraints. Biomaterials 23, 4095–4103. doi:10.1016/S0142-9612(02)00148-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Huiskes, R., Van Driel, W. D., Prendergast, P. J., and Søballe, K. (1997). A biomechanical regulatory model for periprosthetic fibrous-tissue differentiation. J. Mat. Sci. Mat. Med. 8, 785–788. doi:10.1023/a:1018520914512

PubMed Abstract | CrossRef Full Text | Google Scholar

Isaksson, H., van Donkelaar, C. C., Huiskes, R., and Ito, K. (2008). A mechano-regulatory bone-healing model incorporating cell-phenotype specific activity. J. Theor. Biol. 252, 230–246. doi:10.1016/j.jtbi.2008.01.030

PubMed Abstract | CrossRef Full Text | Google Scholar

Lacroix, D., and Prendergast, P. J. (2002). A mechano-regulation model for tissue differentiation during fracture healing: Analysis of gap size and loading. J. Biomechanics 35, 1163–1171. doi:10.1016/S0021-9290(02)00086-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Makowski, P., and Kuś, W. (2016). Optimization of bone scaffold structures using experimental and numerical data. Acta Mech. 227, 139–149. doi:10.1007/s00707-015-1421-4

CrossRef Full Text | Google Scholar

Metz, C., Duda, G. N., and Checa, S. (2020). Towards multi-dynamic mechano-biological optimization of 3D-printed scaffolds to foster bone regeneration. Acta Biomater. 101, 117–127. doi:10.1016/j.actbio.2019.10.029

PubMed Abstract | CrossRef Full Text | Google Scholar

Percoco, G., Uva, A. E., Fiorentino, M., Gattullo, M., Manghisi, V. M., and Boccaccio, A. (2020). Mechanobiological approach to design and optimize bone tissue scaffolds 3D printed with fused deposition modeling: A feasibility study. Materials 13, 648. doi:10.3390/ma13030648

CrossRef Full Text | Google Scholar

Perier-Metz, C., Cipitria, A., Hutmacher, D. W., Duda, G. N., and Checa, S. (2022). An in silico model predicts the impact of scaffold design in large bone defect regeneration. Acta Biomater. 145, 329–341. doi:10.1016/j.actbio.2022.04.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Perier-Metz, C., Duda, G. N., and Checa, S. (2021). Initial mechanical conditions within an optimized bone scaffold do not ensure bone regeneration – An in silico analysis. Biomech. Model. Mechanobiol. 20, 1723–1731. doi:10.1007/s10237-021-01472-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Perier-Metz, C., Duda, G. N., and Checa, S. (2020). Mechano-biological computer model of scaffold-supported bone regeneration: Effect of bone graft and scaffold structure on large bone defect tissue patterning. Front. Bioeng. Biotechnol. 8, 585799. doi:10.3389/fbioe.2020.585799

PubMed Abstract | CrossRef Full Text | Google Scholar

Petersen, A., Princ, A., Korus, G., Ellinghaus, A., Leemhuis, H., Herrera, A., et al. (2018). A biomaterial with a channel-like pore architecture induces endochondral healing of bone defects. Nat. Commun. 9, 4430. doi:10.1038/s41467-018-06504-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Pobloth, A.-M., Checa, S., Razi, H., Petersen, A., Weaver, J. C., Schmidt-Bleek, K., et al. (2018). Mechanobiologically optimized 3D titanium-mesh scaffolds enhance bone regeneration in critical segmental defects in sheep. Sci. Transl. Med. 10, eaam8828. doi:10.1126/scitranslmed.aam8828

PubMed Abstract | CrossRef Full Text | Google Scholar

Poh, P. S. P., Valainis, D., Bhattacharya, K., Griensven, M., and Dondl, P. (2019). Optimization of bone scaffold porosity distributions. Sci. Rep. 9, 9170. doi:10.1038/s41598-019-44872-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Prendergast, P. J., Huiskes, R., and Søballe, K. (1997). Biophysical stimuli on cells during tissue differentiation at implant interfaces. J. Biomechanics 30, 539–548. doi:10.1016/S0021-9290(96)00140-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Reichert, J. C., Cipitria, A., Epari, D. R., Saifzadeh, S., Krishnakanth, P., Berner, A., et al. (2012). A tissue engineering solution for segmental defect regeneration in load-bearing long bones. Sci. Transl. Med. 4, 141ra93. doi:10.1126/scitranslmed.3003720

PubMed Abstract | CrossRef Full Text | Google Scholar

Reznikov, N., Boughton, O. R., Ghouse, S., Weston, A. E., Collinson, L., Blunn, G. W., et al. (2019). Individual response variations in scaffold-guided bone regeneration are determined by independent strain- and injury-induced mechanisms. Biomaterials 194, 183–194. doi:10.1016/j.biomaterials.2018.11.026

PubMed Abstract | CrossRef Full Text | Google Scholar

Schlundt, C., Bucher, C. H., Tsitsilonis, S., Schell, H., Duda, G. N., and Schmidt-Bleek, K. (2018). Clinical and research approaches to treat non-union fracture. Curr. Osteoporos. Rep. 16, 155–168. doi:10.1007/s11914-018-0432-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Sturm, S., Zhou, S., Mai, Y.-W., and Li, Q. (2010). On stiffness of scaffolds for bone tissue engineering - a numerical study. J. Biomech. 43, 1738–1744. doi:10.1016/j.jbiomech.2010.02.020

PubMed Abstract | CrossRef Full Text | Google Scholar

The MathWorks Inc (2020). Global optimization toolbox manual. Natick, Massachusetts, United States. Available at: https://de.mathworks.com/help/gads/(Accessed February 18, 2021).

Google Scholar

Wang, Y., Luo, Z., Zhang, N., and Qin, Q. (2016). Topological shape optimization of multifunctional tissue engineering scaffolds with level set method. Struct. Multidiscipl. Optim. 54, 333–347. doi:10.1007/s00158-016-1409-2

CrossRef Full Text | Google Scholar

Wieding, J., Wolf, A., and Bader, R. (2014). Numerical optimization of open-porous bone scaffold structures to match the elastic properties of human cortical bone. J. Mech. Behav. Biomed. Mater. 37, 56–68. doi:10.1016/j.jmbbm.2014.05.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, C., Fang, J., Entezari, A., Sun, G., V Swain, M., Xu, Y., et al. (2021). A time-dependent mechanobiology-based topology optimization to enhance bone growth in tissue scaffolds. J. Biomechanics 117, 110233. doi:10.1016/j.jbiomech.2021.110233

PubMed Abstract | CrossRef Full Text | Google Scholar

Xiao, D., Yang, Y., Su, X., Wang, D., and Luo, Z. (2012). Topology optimization of microstructure and selective laser melting fabrication for metallic biomaterial scaffolds. Trans. Nonferrous Metals Soc. China 22, 2554–2561. doi:10.1016/S1003-6326(11)61500-8

CrossRef Full Text | Google Scholar