Design and properties of 3D scaffolds for bone tissue engineering

Accepted Manuscript Full length article Design and Properties of 3D Scaffolds for Bone Tissue Engineering S. Gómez, M.D. Vlad, J. López, E. Fernández PII: DOI: Reference:

S1742-7061(16)30309-9 http://dx.doi.org/10.1016/j.actbio.2016.06.032 ACTBIO 4304

To appear in:

Acta Biomaterialia

Received Date: Revised Date: Accepted Date:

9 March 2016 24 May 2016 28 June 2016

Please cite this article as: Gómez, S., Vlad, M.D., López, J., Fernández, E., Design and Properties of 3D Scaffolds for Bone Tissue Engineering, Acta Biomaterialia (2016), doi: http://dx.doi.org/10.1016/j.actbio.2016.06.032

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

3D Voronoi BTE scaffolds

DESIGN AND PROPERTIES OF 3D SCAFFOLDS FOR BONE TISSUE ENGINEERING S. Gómez,1 M.D. Vlad,2 J. López1 and E. Fernández1

1

Research Group of Interacting Surfaces in Bioengineering and Materials Science

(InSup), Technical University of Catalonia (UPC), Avda. Diagonal 647, E-08028Barcelona, Spain. 2

Faculty of Medical Bioengineering, “Gr. T. Popa” University of Medicine and Pharmacy, Str. Kogalniceanu 9-13, 700454-Iasi, Romania.

Address for correspondence and reprints: Enrique Fernández, PhD, Department of Materials Science and Metallurgical Engineering, Technical University of Catalonia (UPC), Avda. Diagonal 647, E-08028Barcelona, Spain; Fax No. +34-93-4016706; Phone No. +34-93-4016287; E-mail: [email protected]

1

3D Voronoi BTE scaffolds Abstract In this study, the Voronoi tessellation method has been used to design novel bone like three dimension (3D) porous scaffolds. The Voronoi method has been processed with computer design software to obtain 3D virtual isotropic porous interconnected models, exactly matching the main histomorphometric indices of trabecular bone (trabecular thickness, trabecular separation, trabecular number, bone volume to total volume ratio, bone surface to bone volume ratio, etc.). These bone like models have been further computed for mechanical (elastic modulus) and fluid mass transport (permeability) properties. The results show that the final properties of the scaffolds can be controlled during their microstructure and histomorphometric initial design stage. It is also shown that final properties can be tuned during the design stage to exactly match those of trabecular natural bone. Moreover, identical total porosity models can be designed with quite different specific bone surface area and thus, this specific microstructural feature can be used to favour cell adhesion, migration and, ultimately, new bone apposition (i.e. osteoconduction). Once the virtual models are fully characterized and optimized, these can be easily 3D printed by additive manufacturing and/or stereolitography technologies.

Keywords: bone tissue engineering, porous scaffold, Voronoi tessellation method, computer-aided tissue engineering, computational fluid dynamics, finite elements analysis

2

3D Voronoi BTE scaffolds 1. Introduction Bone is nowadays one of the most transplanted tissues, with an incidence of nearly 15 million fracture cases per year. Autografts and allografts are available options for bone repair. However, related problems such as pain, infection and immune rejection have favoured the development of artificial bone-like scaffolds. In this case, the state of the art has demonstrated that most common scaffolds lack sufficient mechanical strength and vascularisation [1,2]. In fact, it has been concluded that the optimum bone scaffold should copy natural bone properties at all levels (i.e. mechanical, biological, mass transport and microstructure geometry properties) in order to have similar cells penetration and nutrients diffusion as well as optimum maintenance of biodegradation properties in service [3-14]. This has led Bone Tissue Engineering (BTE) to experiment with several manufacturing technologies having the capacity to reproduce to scale the three-dimensional (3D) microstructure of bone. Among these, Rapid Prototyping (RP) has been revealed as the best technology to fabricate 3D irregular and interconnected porous scaffolds by Additive Manufacturing (AM). This method uses STereoLitography (STL) files, in an optical approach (i.e. the most accurate), and builds the scaffold by printing it layer-by-layer with biocompatible and bioresorbable material. In fact, highresolution stereolitography (µSTL) produces scaffolds with 90% porosity and pores with 20 to 1000 µm of size [15-17]. Nevertheless, to build biomimetic bone scaffolds is not an easy task. However, the design process, that needs to be done before the scaffold’s manufacture, fortunately can be guided by computer assisted design (CAD) software [18]. This has been mainly used to model 3D simple scaffold geometries combining standard primitives (cubes, hexagons, spheres, cylinders, etc.) through union, subtraction and intersection Boolean operations (see, Fig. 1.a). In this way, periodic porous structures can be expanded with 3

3D Voronoi BTE scaffolds the help of a specific array software of commands. Similarly, in order to reduce the operating time of the design, several authors have proposed the use of parametric libraries of regular polyhedrals, known as CASTS (Computer-Aided System for Tissue Scaffolds). However, these elemental scaffolds show poor biomechanical and flow properties [19-23]. Another approach is Implicit Surface (IS) modelling, which uses mathematic trigonometric functions to maximise the surface-volume ratio of the scaffold, this is know to have positive effects on cell migration and tissue ingrowth [24-26]. As an example, Fig. 1.b shows specific plane cuts for several IS-3D models of 85% porosity (made with K3DSurf v.0.6.2 open-software (http://k3dsurf.sourceforge.net)) created from Neovius surface (NS), Schwarz’s Primitive (SP), Schwarz’s Diamond (SD), Schoen’s Gyroid (SG) and Schwarz ‘W’ models, which are the most common to reproduce the porous structure of trabecular bone. As can be observed in Fig. 1.c, the IS-scaffolds of Fig. 1.b show different regular microstructures depending on angular rotation (45º in Fig. 1.c) and, consequently, these models have anisotropic mechanical and fluid flow properties. Additionally, computed tomography (CT) and micro computed tomography (µCT) image analysis have been used as the basis for the design of 3D bio-mimetic scaffolds with irregular

porous structure [27-29]. However,

such complex

microstructure requires expensive facilities, which are often not necessary. Fig. 1.d is an example of how this process works. The figure shows the reconstruction of a µCT-file for the L3 human vertebra [30]. In this article, we present the work dealing with a novel approach, the Voronoi tessellation method, which tries to solve the intrinsic regularities obtained by previous

4

3D Voronoi BTE scaffolds methods. The article shows the suitability of the method when combined with the use of appropriated CAD software. It shows different ways of designing biomimetic bone scaffolds with gradient-controlled interconnected porosity, controlled pore distribution and controlled pore shapes. It also shows how to shape the scaffold to fit real bone defects exactly. Additionally, the paper presents the mechanical and fluid flow properties obtained for several designed scaffolds, over a range of porosities. 2. Materials and methods 2.1 Design and parametric characterization of 3D scaffolds for BTE 3D trabecular bone-like structures were designed using the novel CAD software GrasshopperTM

v.0.9.0076

(http://www.grasshopper3d.com/)

and

the

Voronoi

tessellation method [31]. The design procedure started with the definition of a volumeof-interest (VOI) of 23 mm3 and the inside allocation of certain distribution of points from which to obtain irregular polyhedral unit cells. Two different approaches were followed to define the initial distribution of points. The first one, the most direct, used the histomorphometric bone index “mean trabecular separation” (TbSp) to obtain for certain VOI the mean number (N) of polyhedral cells filling that VOI (i.e. N=[(VOI)/(TbSp)3]) and, consequently, the total number of points that GrasshopperTM distributes in random manner into the VOI [31]. The same software also identifies the spatial coordinates of each of these points, which are necessary for further processing, as it is later explained. The second approach took profit from the µCT images of the L3 human vertebra ESA29-99-L3 (980 slide images with 2048x2048 of pixel resolution and 37 µm of pass resolution)

[30,32].

The

images

were

exported

to

ImageJ©

v.1.47

(http://imagej.nih.gov/) [33] software (see, Fig. 2.a), converted to 8-bit grey scale and, 5

3D Voronoi BTE scaffolds binarized and softened (to minimize noise) with Threshold and Smooth software commands (Fig. 2.b). Then after, the Delaunay-Voronoi command (J. Schindelin and L.P. Chew) for ImageJ© split the images by lines of points having equal distance to the borders of the two nearest particles (Fig. 2.c). Thus, the Voronoi cell that was associated to each particle included all points that were nearer to this particle than any other particle. At this stage, the ImageJ© software command Ultimate Points defined the centre point of each Voronoi cell (Fig. 2.d). To summarize, the process started with 2D images having certain particle structure (µCT slides; Fig. 2.a) and ended with a distribution of single points in 2D (Fig. 2.d), one for each image. The final 3D distribution of points was processed with Rhinoceros3D© software (Robert McNeell & Associates, v.4.0 SR9) summing up consecutive 2D slides at a distance equal to the bone index Tb.Sp. Once the 3D coordinates (x, y and z) of all points are known (see, Fig. 3.a), these are processed by the 3D Voronoi GrasshopperTM software command to obtain the 3D Voronoi cell structure (see, Fig. 3.b). Then, the Explode command is applied to break points, edges, and faces in the 3D Voronoi cells to allow equidistant copy of the faces to create polyhedron cells separated each other an equivalent distance to bone Trabecular Thickness (Tb.Th; see, Fig. 3.c). Then, Boolean operations between the internal and the external contours gave the final 3D porous interconnected structure (see, Fig. 3.d for the external contour and Fig. 3.e for both, the external and internal contours). This structure was further improved with Weaverbird v.0.90.1 (http://www.giuliopiacentino.com/), a plugin for GrasshopperTM software, with the Catmull-Clark smoothing command, which softened the trabecular mesh model (see, Fig. 3.f). It is worth it to mention that the final scaffold’s porosity (related to the bone volume to total volume ratio; BV/TV) can be modified at convenience by changing the 6

3D Voronoi BTE scaffolds Tb.Sp and/or Tb.Th bone indices (see, Fig. 4.a). Similarly, the bone surface to bone volume (BS/BV) ratio can be modified at convenience at constant porosity by controlling the number of the nucleating points (Fig. 4.b). These new design-related possibilities should find profit to understand ageing of bone tissue properties [34]. In fact, the explained methodology is so powerful that it allows the definition of porous models fitting real bone defects (Fig. 5) and even the definition of porous models with different porosity gradients along different directions so as to match real bone structures (Fig. 6). 2.2 Mechanical characterization of 3D scaffolds for BTE 21 scaffold models were generated, as explained before, using the first approach, i.e. random distributing 10, 15 and 20 points inside the VOI, and designing for the histomorphometric bone indices: 10 < BV/TV(%) < 65; 100 < Tb.Th(µm) < 400; and, 25 < Tb.Sp(µm) < 225. The obtained 3D models were converted to closed solid polysurfaces (mean value of 20000 surfaces obtained from closed polygon meshes; see, Figs. 7.a-b) with the Meshtonurbs command of Rhinoceros3D© and exported to Comsol Multiphysics© v.4.2.a software as IGES (Initial Graphics Exchange Specification) files, for their accurate mechanical computation (Hewlett-Packard© Workstation: Intel© Core i5, 96 GB; 350000 tetrahedral mesh elements on average; maximum time for meshing: 10 min; and minimum time for analysis: 5 min). The homogeneous, isotropic and linear elastic material selected for computation was poly(D,L-lactide), one of the most frequently used bone tissue engineering materials for 3D printing (3DP). The properties selected for this material were: Young’s modulus, Es = 3.3 GPa; Poisson’s ratio, ν = 0.3; and, density, ρ = 1.3 g/cm3 [6,35,36].

7

3D Voronoi BTE scaffolds Compression experiments were simulated (see, Fig. 7.c) under the following boundary conditions: displacement for the top-surface (maximum engineering strain, εmax = 1 %); fixation for the bottom-surface; and, symmetry for the rest of the surfaces [35]. The simulations reported for the maximum (σmax) and the average true reaction stress at the bottom surface (σav) as well as for the average von Mises stress (σVM), displacement field (df) and strain energy density (Estrain). Additionally, the Young’s modulus of the porous material Ep was calculated as the ratio between the σav and the normalized df (or axial engineering strain) at a value of 0.005 (i.e. 0,5%). 2.3 Fluid flow characterization of 3D scaffolds for BTE The fluid flow computation domain was obtained in Rhinoceros3D© after Boolean subtraction of the scaffold geometry (see, Fig. 8.a) from the volume-of-interest. These domains (see, Figs. 8.b-c) were exported, as explained before, in IGES format to Comsol Multiphysics©. Computation Fluid Dynamic (CFD) analysis of the scaffolds was approximated by laminar and stationary Navier-Stokes model. The simulations were done along the x, y and z axes to account for anisotropy (see, Figs. 8.d-f). The boundary conditions defined an inlet-flow side (inlet velocity, vi = 1 mm/s) and an output opposite-flow side (zero output pressure), with cube lateral faces defined symmetric and inner faces as non-slip [37]. Dulbecco's modified Eagle's medium (DMEM) was the fluid material (dynamic viscosity, µd = 1.45 mPa.s; density, ρ = 1 g/cm3 [38]) selected for computation. The convergence analysis concluded that 1.2 million tetrahedral elements on average were required for accurate computation of results. The results reported the average pressure drop between the inlet and the outlet (∆Pi-o(Pa)) faces of area A, and the average velocity (vf(m/s)), shear rate (γ(s-1)) and

8

3D Voronoi BTE scaffolds vorticity (VOR(s-1)) of the fluid across the length L of the model (L = (VOI)1/3). These values allow calculation of scaffold’s intrinsic permeability (K(m2)) according to Darcy's law (Eq. 1) [39-43], where vD(m/s) is the Darcy’s velocity, defined as the volumetric flow rate Q(m3/s) per unit area A(m2) of the sample, and the other parameters as defined previously; in our case, vD = vf. ௅

‫ݒ = ܭ‬஽ ∙ ߤௗ ∙ ቀ∆௉

೔ష೚

ቁ

[Eq. 1]

3. Results 3.1 Design and parametric characterization of 3D scaffolds for BTE Table 1, which also contains the main abbreviations to easily allow the reader to follow the article, summarizes the main histomorphometric bone indices obtained for certain Voronoi cellular structures having similar graded porosity but different initial number of nucleating points (NNP). It is observed that for a fixed value of NNP a survey of different porous structures can be obtained in such a way that when BV/TV increases, BS/TV and Tb.Th also increase but Tb.Sp and Tb.N decrease (see also Figs. 9.a. and 9.b). It is also observed that for increasing values of NNP but similar values of BV/TV, Tb.Th and Tb.Sp decrease but BS/TV and Tb.N increase, as expected. In addition, Fig. 10 summarizes the functionality of the BS/BV index versus the porosity (i.e. 1-BV/TV) for several Voronoi scaffold models and for the most common IS-models (see also Figs. 1.b. and 1.c). 3.2 Mechanical characterization of 3D scaffolds for BTE Figs. 11.a-b show, for some Voronoi scaffold models (NNP = 10, 15 and 20), the evolution of the normalized Young’s modulus (Ep/Es), i.e. the ratio between the elastic

9

3D Voronoi BTE scaffolds modulus of both, the Voronoi porous structure (Ep) and the polymeric reference material (Es), versus the histomorphometric indices Tb.Th and BV/TV, respectively. The main observation is that the average elastic modulus () increases linearly, as can be seen in the correlation (see continuous black line and fitting equations in the figures), with both indices. 3.3 Fluid flow characterization of 3D scaffolds for BTE Fig. 12 shows the relation between the average (x, y and z) computed permeability and the total porosity (i.e. Φ=1-BV/TV) both independently measured for different Voronoi scaffold models (i.e. NNP= 10, 15 and 20). Fig. 13 shows the good correlation observed between the specific surface area and porosity. 4. Discussion In this study, the Voronoi tessellation method has been used to design novel bone like three dimension (3D) porous scaffolds by processing with computer design software to obtain 3D virtual isotropic porous interconnected models, exactly matching (see Table 1) the main histomorphometric indices of trabecular bone [32]. In this sense, several observations can be made from our study. Consequently, Fig. 10 shows that if BS/BV is constant (and so, Tb.Th (see abbreviations in Table 1)), higher porosity Voronoi models can be obtained for lower values of the number of nucleating points, NNP, as expected. Moreover, it is also observed that by adjusting the NNP an appropriate Voronoi scaffold model can be designed to adjust the properties of its equivalent Implicit Surface model. For example, a Voronoi scaffold model with 10 NNP gives nearly the same results than both, Cylinder Grid and Schwarz Primitive models. Similarly, the studied 100 NNP or the 125 NNP Voronoi scaffold models are in between the Schwarz ‘W’ and the Schoen Gyroid or de Neovius’ Surface and Schwarz

10

3D Voronoi BTE scaffolds Diamond models, respectively. In any case, if porosity is set to certain fixed value the model that gives the higher value of BS/BV index is the Neovius’ Surface model while the one that gives the lower value is the Schwarz ‘W’. However, a Voronoi scaffold model can always be designed to equal those Implicit Surface models. In this sense, the proposed Voronoi methodology indicates that by fixing Tb.Th, Tb.Sp and NNP suitable 3D bone-like porous scaffolds can be created to mimic structural features of natural bones. Moreover, as the Tb.Th and BV/TV indices are also interrelated (see Fig. 9), the general conclusion is again the same as before, i.e. that suitable mechanical and structural bone-like isotropic porous scaffold can be designed by an appropriate selection of the indices. This means that for certain fixed BV/TV value of interest, the Voronoi scaffold with the highest BS/TV ratio (see, Fig. 9.b; i.e. higher NNP value) is that with also the lowest elastic modulus of the Voronoi porous scaffold Ep (see, Fig. 11.b), Tb.Th (see, Figs. 9.a and 11.a) and Tb.Sp (see Fig. 9.a). Consequently, the scaffold with the highest stiffness is that generated with the lower number of nucleating points value NNP and the highest BV/TV ratio (see, Fig. 11.b). In this case, the 3D scaffold has the highest Tb.Th (see, Fig. 9.a) and the lowest BS/TV ratio (see, Fig. 9.b). Additionally, the results also show that the scaffold Young’s modulus is not significantly different for BV/TV lower than 15% (i.e. porosity higher than 85%). However, in this range, the scaffolds show an increasing surface area (i.e. BS/TV increases) when increasing the number of nucleating points, NNP. For this reason, it is expected that bone tissue engineering scaffolds designed as above should show increasing rates of cellular adhesion. The isotropic performance of the designed scaffolds was also confirmed. In fact, compressive studies performed along the three-symmetry axis showed that isotropy 11

3D Voronoi BTE scaffolds increased with NNP and that no relevant differences were observed for NNPs higher than 20. On the other hand, our compressive results (not shown) were in agreement with those measured for trabecular bone (i.e. 20 to 500 MPa) for similar porosity values [35,36,44-49]. The average permeability showed (Fig. 12), through its standard deviation, that models could be considered isotropic. Each scaffold model in Fig. 10 has been fitted to the classical power law model (K∼Φn) to show the main trend (correlation coefficients, 0.86
12

3D Voronoi BTE scaffolds The results obtained in this study show that both, the mechanical and the mass transport properties of Voronoi implants can be properly balanced; while mechanical strength increases at the expense of porosity reduction, permeability increases (for interconnected porosity) when specific bone surface area (and consequently, the number of nucleating points) is decreased. In this sense, the study shows that Voronoi scaffold models having porosity of 75-85% will also have elastic modulus of ∼300-500 MPa, trabecular thickness of ∼125-225 µm, trabecular separation of ∼125-200 µm, bone specific surface area of 15-30 mm-1 and permeability of ∼1.5-3.5x10-8 m2; all these values in the good order to favour diffusion of essential nutrients and oxygen for cell survivability and ultimately osteogenesis [35,36,44-46,50,51]. In addition, the scaffolds’ isotropic properties facilitate their insertion into bone defects regardless of its orientation. Thus, the final shape of an implant scaffold is obtained by Boolean difference operation between the implant and the anatomical shape of the defect to be filed. Finally, the proposed methodology can also create scaffolds with varying but connected porosity in different volumes of interest. Thus, by using medical reverse engineering technology it is possible to create, for example, complex human trabecular bone scaffold models from three-dimensional reconstruction of tomography images. Then, the STL files can be used to print the 3D scaffold with biocompatible o bioinert material by additive manufacturing techniques as stereolithography. These scaffolds can be further treated with specific cells and nutrients for specific bone tissue engineering applications. 5. Conclusion The present study shows that the Voronoi tessellation method is a versatile mathematical method that linked to both, appropriate computer assisted design software and additive manufacturing technology, can favour the ultimate fabrication of mimic

13

3D Voronoi BTE scaffolds bone like implantable structures, copying natural bone properties at all levels (i.e. microstructure, mechanical, mass transport and biological properties), with optimum cell penetration, nutrients diffusion and osteoconduction properties. The results have shown that the computed mechanical and fluid properties of Voronoi scaffold models directly depend on the microstructural features of the porous structure, being total porosity and bone surface area the main factors to control, these depending also on microstructural fundamental bone indices like trabecular thickness, trabecular separation and trabecular number, all of these indices easily controlled at the very beginning of the Voronoi design methodology process that has been proposed. Future research should focus in the manufacture of the 3D Voronoi scaffolds and their in vitro and in vivo behaviour.

14

3D Voronoi BTE scaffolds Acknowledgements The authors thank public funding received for this work through project MAT201019431 (Ministerio de Ciencia e Innovación, Spain). A special thanks to Ms. Mandy Farrell for final editing of the manuscript in English.

References [1] R.J. O'Keefe, J. Mao, Bone tissue engineering and regeneration: from discovery to the clinic - an overview, Tissue Eng Part B Rev 17 (2011) 389-392. https://doi.org/10.1089/ten.TEB.2011.0475. [2] Y. Liu, J. Lim, S.H. Teoh, Review: development of clinically relevant scaffolds for vascularised

bone

tissue engineering,

Biotechnol Adv 31

(2013) 688-705.

https://doi.org/10.1016/j.biotechadv.2012.10.003. [3] S. Bose, M. Roy, A. Bandyopadhyay, Recent advances in bone tissue engineering scaffolds,

Trends

Biotechnol

30

(2012)

546-554.

https://doi.org/10.1016/j.tibtech.2012.07.005. [4] D.W. Hutmacher, J.T. Schantz, C.X. Lam, K.C. Tan, T.C. Lim, State of the art and future directions of scaffold-based bone engineering from a biomaterials perspective, J Tissue Eng Regen Med 1 (2007) 245-260. https://doi.org/10.1002/term.24. [5] B.M. Willie, A. Petersen, K. Schmidt-Bleek, A. Cipitria, M. Mehta, P. Strube, J. Lienau, B. Wildemann, P. Fratzl, G. Duda, Designing biomimetic scaffolds for bone regeneration: why aim for a copy of mature tissue properties if nature uses a different approach?, Soft Matter 6 (2010) 4976-4987. https://doi.org/10.1039/C0SM00262C. [6] D.W. Hutmacher, Scaffolds in tissue engineering bone and cartilage, Biomaterials 21 (2000) 2529-2543. https://doi.org/10.1016/S0142-9612(00)00121-6.

15

3D Voronoi BTE scaffolds [7] S.J. Hollister, Porous scaffold design for tissue engineering, Nat Mater 4 (2005) 518-524. https://doi.org/10.1038/nmat1683. [8] K. Rezwan, Q.Z. Chen, J.J. Blaker, A.R. Boccaccini, Biodegradable and bioactive porous polymer/inorganic composite scaffolds for bone tissue engineering, Biomaterials 27 (2006) 3413-3431. https://doi.org/10.1016/j.biomaterials.2006.01.039. [9] Y. Khan, M.J. Yaszemski, A.G. Mikos, C.T. Laurencin, Tissue engineering of bone: material and matrix considerations, J Bone Joint Surg Am 90 (2008) 36-42. https://doi.org/10.2106/JBJS.G.01260. [10] K.F. Leong, C.M. Cheah, C.K. Chua, Solid freeform fabrication of threedimensional scaffolds for engineering replacement tissues and organs, Biomaterials 24 (2003) 2363-2378. https://doi.org/10.1016/S0142-9612(03)00030-9. [11] P. Lichte, H.C. Pape, T. Pufe, P. Kobbe, H. Fischer, Scaffolds for bone healing: concepts,

materials

and

evidence,

Injury

42

(2011)

569-573.

https://doi.org/10.1016/j.injury.2011.03.033. [12] X. Liu, P.X. Ma, Polymeric scaffolds for bone tissue engineering, Ann Biomed Eng 2004 (32) 477-486. Doi: 10.1023/B:ABME.0000017544.36001.8e [13] V. Karageorgiou, D. Kaplan, Porosity of 3D biomaterial scaffolds and osteogenesis,

Biomaterials

26

(2005)

5474-5491.

https://doi.org/10.1016/j.biomaterials.2005.02.002. [14] L.M. Mathieu, T.L. Mueller, P.E. Bourban, D.P. Pioletti, R. Müller, J.E. Månson, Architecture and properties of anisotropic polymer composite scaffolds for bone tissue engineering,

Biomaterials

27

(2006)

905-916.

https://doi.org/10.1016/j.biomaterials.2005.07.015.

16

3D Voronoi BTE scaffolds [15] F.W. Melchels, M.N. Domingos, T.J. Klein, J. Malda, P.J. Bartolo, D.W. Hutmacher, Additive manufacturing of tissues and organs, Prog Polym Sci 37 (2012) 1079-1104. https://doi.org/10.1016/j.progpolymsci.2011.11.007. [16] F.W. Melchels, J. Feijen, D.W. Grijpma, A review on stereolithography and its applications in biomedical engineering. Biomaterials 31 (2010) 6121-6130. Doi: https://doi.org/10.1016/j.biomaterials.2010.04.050. [17] K. Kim, A. Yeatts, D. Dean, J.P. Fisher, Stereolithographic bone scaffold design parameters: osteogenic differentiation and signal expression, Tissue Eng Part B Rev 16 (2010) 523-539. https://doi.org/10.1089/ten.teb.2010.0171. [18] Computer-Aided Tissue Engineering. Michael A.K. Liebschner (Ed.). Humana Press, Springer Science+Business Media, LLC, 233 Spring Street, New York, NY 10013, USA (2012). ISBN 978-1-61779-763-7. https://doi.org/10.1007/978-1-61779764-4. [19] S.M. Giannitelli, D. Accoto, M. Trombetta, A. Rainer, Current trends in the design of scaffolds for computer-aided tissue engineering, Acta Biomater 10 (2014) 580-594. https://doi.org/10.1016/j.actbio.2013.10.024. [20] D.W. Hutmacher, M. Sittinger, M.V. Risbud, Scaffold-based tissue engineering: rationale for computer-aided design and solid free-form fabrication systems, Trends Biotechnol 22 (2004) 354-362. https://doi.org/10.1016/j.tibtech.2004.05.005. [21] W. Sun, B. Starly, A. Darling, C. Gomez, Computer-aided tissue engineering: application to biomimetic modelling and design of tissue scaffolds, Biotechnol Appl Biochem 39 (2004) 49-58. https://doi.org/10.1042/BA20030109.

17

3D Voronoi BTE scaffolds [22] C.M. Cheah, C.K. Chua, K.F. Leong, S.W. Chua, Development of a tissue engineering scaffold structure library for rapid prototyping. Part 1: Investigation and classification,

Int

J

Adv

Manuf

Technol

21

(2003)

291-301.

https://doi.org/10.1007/s001700300034. [23] C.M. Cheah, C.K. Chua, K.F. Leong, S.W. Chua, Development of a tissue engineering scaffold structure library for rapid prototyping. Part 2: Parametric library and

assembly program,

Int

J

Adv Manuf

Technol 21

(2003)

302-312.

https://doi.org/10.1007/s001700300035. [24] S.C. Kapfer, S.T. Hyde, K. Mecke, C.H. Arns, G.E. Schröder-Turk, Minimal surface scaffold designs for tissue engineering, Biomaterials 32 (2011) 6875-6882. https://doi.org/10.1016/j.biomaterials.2011.06.012. [25] D.J. Yoo, Computer-aided porous scaffold design for tissue engineering using triply periodic minimal surfaces, International Journal of Precision Engineering and Manufacturing 12 (2011) 61-71. https://doi.org/10.1007/s12541-011-0008-9. [26] W. Sun, P. Lal, Recent development on computer aided tissue engineering - a review,

Comput

Methods

Programs

Biomed

67

(2002)

85-103.

https://doi.org/10.1016/S0169-2607(01)00116-X. [27] Z. Chen, Z. Su, S. Ma, X. Wu, Z. Luo, Biomimetic modeling and three-dimension reconstruction of the artificial bone, Comput Methods Programs Biomed 88 (2007) 123130. https://doi.org/10.1016/j.cmpb.2007.08.001. [28] W. Sun, B. Starly, J. Nam, A. Darling, Bio-CAD modeling and its applications in computer-aided tissue engineering, Comput Aided Des 37 (2005) 1097-1114. https://doi.org/10.1016/j.cad.2005.02.002.

18

3D Voronoi BTE scaffolds [29] X.Y. Kou, S.T. Tan, A simple and effective geometric representation for irregular porous structure

modeling,

Comput Aided Des 42

(2010) 930-941. Doi:

https://doi.org/10.1016/j.cad.2010.06.006. [30] G. Beller, M. Burkhart, D. Felsenberg, W. Gowin, H.C. Hege, B. Koller, S. Prohaska, P.I. Saparin, J.S. Thomsen. Vertebral Body Data Set ESA29-99-L3, 2005. [31] A. Okabe, B. Boots, K. Sugihara, S.N. Chiu, Spatial Tessellations: Concepts and Applications of Voronoi Diagrams, 2nd Edition, John Wiley & Sons Ltd, Baffins Lane, England, 2000. ISBN 0-471-98635-6. [32] S. Gómez, M.D. Vlad, J. López, M. Navarro, E. Fernández, Characterization and three-dimensional reconstruction of synthetic bone model foams, Mater Sci Eng C Mater Biol Appl 33 (2013) 3329-3335. https://doi.org/10.1016/j.msec.2013.04.013. [33] W.S. Rasband, ImageJ, U. S. National Institutes of Health, Bethesda, Maryland, USA, version 1.47p, http://imagej.nih.gov/ij/, 1997-2012. [34] A.M. Parfitt, C.H. Mathews, A.R. Villanueva, M. Kleerekoper, B. Frame, D.S. Rao, Relationships between surface, volume, and thickness of iliac trabecular bone in aging and in osteoporosis. Implications for the microanatomic and cellular mechanisms of bone loss, J Clin Invest 72 (1983) 1396-1409. https://doi.org/10.1172/JCI111096. [35] M.R.A. Kadi, A. Syahrom, A. Öchsner, Finite element analysis of idealised unit cell trabecular structure based on morphological indices of cancellous bone, Medical & Biological

Engineering

&

Computing

48

(2010)

497-505.

https://doi.org/10.1007/s11517-010-0593-2. [36] L.J. Gibson, Biomechanics of cellular solids, J Biomech 38 (2005) 377-399. https://doi.org/10.1016/j.jbiomech.2004.09.027.

19

3D Voronoi BTE scaffolds [37] A. Syahrom, M.R.A. Kadir, J. Abdullah, A. Öchsner, Permeability studies of artificial and natural cancellous bone structures, Med Eng Phys 35 (2013) 792-799. https://doi.org/10.1016/j.medengphy.2012.08.011. [38] A.L. Olivares, E. Marsal, J.A. Planell, D. Lacroix, Finite element study of scaffold architecture design and culture conditions for tissue engineering, Biomaterials 30 (2009) 6142-6149. https://doi.org/10.1016/j.biomaterials.2009.07.041. [39] R.P. Widmer, S.J. Ferguson, A comparison and verification of computational methods to determine the permeability of vertebral trabecular bone, Proc Inst Mech Eng H 227 (2013) 617-628. https://doi.org/10.1177/0954411912462814. [40] E.A. Nauman, K.E. Fong, T.M. Keaveny, Dependence of intertrabecular permeability on flow direction and anatomic site, Ann Biomed Eng 27 (1999) 517-524. https://doi.org/10.1114/1.195. [41] J.F. Despois, A. Mortensen, Permeability of open-pore microcellular materials, Acta Mater 53 (2005) 1381-1388. https://doi.org/10.1016/j.actamat.2004.11.031. [42] R. Singh, P.D. Lee, T.C. Lindley, R.J. Dashwood, E. Ferrie, T. Imwinkelried, Characterization of the structure and permeability of titanium foams for spinal fusion devices, Acta Biomater 5 (2009) 477-487. https://doi.org/10.1016/j.actbio.2008.06.014. [43] M.J. Grimm, J.L. Williams, Measurements of permeability in human calcaneal trabecular bone, J Biomech 30 (1997) 743-745. https://doi.org/10.1016/S00219290(97)00016-X. [44] D. Ulrich, B. Rietbergen, A. Laib, P. R̈ uegsegger, The ability of three-dimensional structural indices to reflect mechanical aspects of trabecular bone, Bone 25 (1999) 5560. https://doi.org/10.1016/S8756-3282(99)00098-8.

20

3D Voronoi BTE scaffolds [45] S. Majumdar, M. Kothari, P. Augat, D.C. Newitt, T.M. Link, J.C. Lin, T. Lang, Y. Lu, H.K. Genant, High-resolution magnetic resonance imaging: three-dimensional trabecular bone architecture and biomechanical properties, Bone 22 (1998) 445-454. https://doi.org/10.1016/S8756-3282(98)00030-1. [46] D.W. Hutmacher, Scaffolds in tissue engineering bone and cartilage, Biomaterials 21 (2000) 2529-2543. https://doi.org/10.1016/S0142-9612(00)00121-6. [47] R.W. Goulet, S.A. Goldstein, M.J. Ciarelli, J.L. Kuhn, M.B. Brown, L.A. Feldkamp, The relationship between the structural and orthogonal compressive properties

of

trabecular

bone,

J

Biomech

27

(1994)

375-389.

https://doi.org/10.1016/0021-9290(94)90014-0. [48] T.A. Burgers, J. Mason, G. Niebur, H.L. Ploeg, Compressive properties of trabecular

bone

in

the

distal

femur,

J

Biomech

41

(2008)

1077-1085.

https://doi.org/10.1016/j.jbiomech.2007.11.018. [49] S.A. Goldstein, The mechanical properties of trabecular bone: dependence on anatomic

location

and

function,

J

Biomech

20

(1987)

1055-1061.

https://doi.org/10.1016/0021-9290(87)90023-6. [50] A.E. Johnson, T.S. Keller, Mechanical properties of open-cell foam synthetic thoracic

vertebrae,

J

Mater

Sci

Mater

Med

19

(2008)

1317-1323.

https://doi.org/10.1007/s10856-007-3158-7. [51] J.H. Kinney, A.J.C. Ladd, The relationship between three-dimensional connectivity and the elastic properties of trabecular bone, J Bone Miner Res 13 (1998) 839-845. https://doi.org/10.1359/jbmr.1998.13.5.839.

21

3D Voronoi BTE scaffolds

Figure Legends Figure 1. a) Different scaffold porous units obtained from Boolean operations of standard primitives (cubes, hexagons, spheres, cylinders, etc.). b) Specific internal plane cuts of 3D implicit surface regular 85% porous models; from left to right: Neovius surface (NS), Schwarz primitive (SP), Schwarz diamond (SD), Schoen Gyroid (SG) and Schwarz ‘W’. c) Similar as above but after 45º of plane rotation to show the intrinsic anisotropy of these IS models. d) From left to right: selection of the region of interest (ROI) of consecutive µCT images of human vertebra; binarization; smoothing; and 3D reconstruction of the ROI. Figure 2. Process to obtain 2D Voronoi cell structures from µCT images of L3 human vertebra. a) µCT images are exported to ImageJ© software and converted to 8-bit grey scale. b) Then after, the images are binarized and softened with the Threshold and Smooth ImageJ© software commands. c) Then after, the images are split with the Delaunay-Voronoi command for ImageJ© software. d) Finally, the Voronoi cell structure is obtained with the Ultimate Points ImageJ© software command. Figure 3. Process to obtain 3D Voronoi cell regular porous structures. a) 2D Voronoi point coordinates (x,y) are processed at equal z distances to obtain the 3D coordinates (x,y,z) of all points. b) Then after, the coordinates are processed to obtain the 3D Voronoi cell structure. c) Then after, each plane surface is self-copy and translated a fix distance to create an open internal connected volume. d) and e) Then after, the porous external (c) and internal (d) contours are obtained by Boolean operations. f) Finally, a smoothing process is done to soften the trabecular mesh model.

22

3D Voronoi BTE scaffolds Figure 4. a) The Voronoi scaffolds’ porosity is controlled for a fixed number of nucleating points (NNP) by changing both trabecular separation Tb.Sp and thickness Tb.Th indices. b) Similarly, the scaffolds’ specific surface area is changed at constant porosity Φ by adjusting the NNPs. Figure 5. Steps followed to model a 3D porous scaffold fitting both an external and internal predefined shape. a) Selection and definition of the implantable volume zone. b) Virtual 3D model of the implant. c) Modification of the inner part of the implant to be porous-like trabecular bone. Figure 6. a) Porous scaffold designed with a gradient of open and interconnected porosity. The central part of the implant has higher porosity than the external parts. This design is more self-similar to real bone microstructures. b) 3D microstructure of human vertebra ESA29-99-L3 [30] where gradients of porosity are evident (piece of 18x8x8 mm3). Figure 7. Finite element procedure used for the mechanical characterization of the scaffolds. a) Plane projection of the 3D scaffold. The computation conditions where: bottom and top faces fixed and displaced in the y-direction, respectively; maximum strain for the top face of 0.5 %; lateral faces all symmetric. b) Detail showing the tetrahedral fine element mesh to obtain optimum convergence computation. c) Example of a contour colour plot obtained for the von Mises stress results. Figure 8. Steps followed for the fluid dynamic characterization of the scaffolds. a) A porous scaffold is selected. b) The fluid computation domain is obtained by Boolean subtraction operation. c) An appropriate convergent meshing model is obtained. d-e) Computation is performed along the three axes to account for anisotropy (from left to right: x, y and z axis).

23

3D Voronoi BTE scaffolds Figure 9. a) Trabecular thickness (Tb.Th) and trabecular separation (Tb.Sp) indices as a function of the bone volume to total volume (BV/TV) ratio for three Voronoi scaffold models with 10, 15 and 20 nucleating points. b) Bone surface to total volume (BS/TV) ratio and trabecular number (Tb.N) indices as a function of BV/TV ratio for the same above models. (Note: Consider Ref. [32] to compare these values to those of natural bone) Figure 10. Bone surface to bone volume (BS/BV) ratio index as a function of porosity for several implicit surface models as compared to some Voronoi scaffold models. (Note: Consider Ref. [32] to compare these values to those of natural bone) Figure 11. Normalized Young’s modulus results obtained for several Voronoi scaffold models with different number of nucleating points as a function of: a) the trabecular thickness (Tb.Th) index, and b) the bone volume to total volume (BV/TV) ratio index. (Note: Consider Refs. [44-49] to compare these values to those of natural bone) Figure 12. Permeability versus porosity for some Voronoi scaffold models. (Note: Consider Refs. [37,39-43] to compare these values to those of natural bone) Figure 13. Specific surface area or bone surface to bone volume (BS/BV) ratio index versus porosity for the studied Voronoi scaffold models. (Note: Consider Ref. [32] to compare these values to those of natural bone)

24

3D Voronoi BTE scaffolds

Figure 1

(a)

NS

SP

SD (b)

SG

SW

NS-45º

SP-45º

SD-45º (c)

SG-45º

SW-45º

µCT image stack

ROI - Binarized

ROI - Smoothed

(d)

3D reconstruction

3D Voronoi BTE scaffolds

Figure 2

(a)

(b)

(c)

(d)

3D Voronoi BTE scaffolds

Figure 3

(a)

(b)

(c)

(d)

(e)

(f)

3D Voronoi BTE scaffolds

Figure 4

(i) NNP = 5; Φ(%) = 89

(ii) NNP = 5; Φ(%) = 69 (a)

(iii) NNP = 5; Φ(%) = 49

(i) NNP = 5; Φ(%) = 89

(ii) NNP = 25; Φ(%) = 89 (b)

(iii) NNP = 50; Φ(%) = 89

3D Voronoi BTE scaffolds

Figure 5

(a)

(b)

(c)

3D Voronoi BTE scaffolds

Figure 6

Φi

Φj > Φi

(a)

Φi

(b)

3D Voronoi BTE scaffolds

Figure 7

(a)

(b)

(c)

3D Voronoi BTE scaffolds

Figure 8

(a)

(b)

(c)

(d)

(e)

(f)

3D Voronoi BTE scaffolds

Figure 9 425 375 325

Tb.Th 10 Tb.Th 15 Tb.Th 20 Tb.Sp (µm)

Tb.Th (µm)

275 225 175

Tb.Sp 10 Tb.Sp 15 Tb.Sp 20

125 75 25 10

20

30

40

50

60

70

BV/TV (%)

(a)

8

7

BS/TV 10 BS/TV 15 BS/TV 20 Tb.N (mm )

-1

-1

BS/TV (mm )

6

5

4

Tb.N 10 Tb.N 15 Tb.N 20

3

2 10

20

30

40

50

60

70

BV/TV (%)

(b)

33

3D Voronoi BTE scaffolds

Figure 10

65

Cylinder Grid Voronoi 010 Schwarts Primitive

45

Schwarts W Voronoi 100 Schoen Gyroid

-1

BS/BV (mm )

55

35

Nervius' Surface Voronoi 125 Schwarts Diamond

25

15

5 30

40

50

60

70

80

90

100

Porosity (%)

34

3D Voronoi BTE scaffolds

Figure 11

0,5

Normalized Young's Modulus (Ep/Es)

= -0,14322 + 0,00147*(Tb.Th) 2

R = 0,957 +/- 0.026 (p < 0.0001)

0,4

0,3

0,2

Voronoi 10 Voronoi 15 Voronoi 20

0,1

0,0 100

150

200

250

300

350

400

Tb.Th (µm)

(a)

0,5

Normalized Young's Modulus (Ep/Es)

= -0,03912 + 0,00714*(BV/TV) 2

R = 0,944 +/- 0.030 (p < 0.0001)

0,4

0,3

0,2

Voronoi 10 Voronoi 15 Voronoi 20

0,1

0,0 0

10

20

30

40

50

60

70

BV/TV (%)

(b)

35

3D Voronoi BTE scaffolds

Figure 12

5,0 4,5

3,5

-8

2

Permeability, K (x10 m )

4,0

Voronoi 10 Voronoi 15 Voronoi 20

3,0 2,5 2,0 1,5 1,0 0,5 0,0 30

40

50

60

70

80

90

100

Porosity, Φ (%)

36

3D Voronoi BTE scaffolds

Figure 13

35

-1

Specific Surface Area, BS/BV (mm )

40

Voronoi 10 Voronoi 15 Voronoi 20

30

25

20

15

10

5 30

40

50

60

70

Porosity, Φ (%)

80

90

100

3D Voronoi BTE scaffolds

Table 1 Table 1. Histomorphometric bone indices of Voronoi cellular structures with similar graded porosity but different number of nucleating points. (Note: Use Ref. 32 to compare these values with those of human vertebral bone) NNPa

BV/TVb (%)

BS/BVc (mm-1)

BS/TVd (mm-1)

Tb.The,h (µm)

Tb.Nf,h (mm-1)

Tb.Spg,h (µm)

12,22 28,62 3,50 140 2,82 215 18,80 21,76 4,09 184 2,66 192 29,66 16,77 4,97 239 2,58 150 10 35,00 14,65 5,13 273 2,44 136 41,28 13,25 5,47 302 2,40 115 53,73 11,11 5,97 360 2,30 75 60,16 10,14 6,10 394 2,22 56 12,89 29,91 3,86 134 3,03 196 18,50 25,22 4,67 159 3,06 168 29,13 18,67 5,44 214 2,84 138 15 34,92 16,54 5,78 242 2,76 121 41,31 15,03 6,21 266 2,73 101 53,43 12,51 6,68 320 2,58 68 61,05 11,23 6,86 356 2,48 48 12,18 35,27 4,30 113 3,47 175 18,32 27,29 5,00 147 3,30 157 29,30 20,17 5,91 198 3,08 126 20 35,06 17,86 6,26 224 2,98 111 40,49 16,69 6,76 240 3,00 94 53,63 13,38 7,18 299 2,76 63 61,51 12,13 7,46 330 2,68 43 a NNP; Number of Nucleating Points. bBV/TV; Bone Volume to Total Volume ratio. cBS/BV; Bone Surface to Bone Volume ratio. dBS/TV=(BS/BV)x(BV/TV); Bone Surface to Total Volume ratio. e Tb.Th=4/(BS/BV); Trabecular Thickness or Trabecular Diameter. fTb.N=((4/π)x(BV/TV))1/2/(Tb.Th); Trabecular Number. gTb.Sp=Tb.N-1 -Tb.Th; Trabecular Separation. h Cylindrical rod model.

38

i

j > i

i

5,0

Cylinder Grid Voronoi 010 Schwarts Primitive

4,5 4,0

35

= -0,03912 + 0,00714*(BV/TV)

3,5

-8

2

Permeability, K (x10 m )

45

Schwarts W Voronoi 100 Schoen Gyroid

-1

BS/BV (mm )

55

0,5

Voronoi 10 Voronoi 15 Voronoi 20

Normalized Young's Modulus (E p/Es)

65

Nervius' Surface Voronoi 125 Schwarts Diamond

25

3,0 2,5 2,0 1,5 1,0

15 0,5

5 30

40

50

60

70

Porosity (%)

80

90

100

0,0 30

40

50

60

70

Porosity,  (%)

80

90

100

2

R = 0,944 +/- 0.030 (p < 0.0001)

0,4

0,3

0,2

Voronoi 10 Voronoi 15 Voronoi 20

0,1

0,0 0

10

20

30

40

BV/TV (%)

50

60

70

Statement of Significance The significance of this article goes far beyond the specific objectives on which it is focussed. In fact, it shows, in a guided way, the entire novel process that can be followed to design graded porous implants, whatever its external shape and geometry, but internally tuned to the exact histomorphometric indices needed to match natural human tissues microstructures and, consequently, their mechanical and fluid properties, among others. The significance is even more relevant nowadays thanks to the available new computing and design software that is easily linked to the 3D printing new technologies. It is this transversality, at the frontier of different disciplines, the main characteristic that gives this article a high scientific impact and interest to a broaden audience.