Height-to-Diameter Ratio and Porosity Strongly Influence Bulk Compressive Mechanical Properties of 3D-Printed Polymer Scaffolds

Although the architectural design parameters of 3D-printed polymer-based scaffolds—porosity, height-to-diameter (H/D) ratio and pore size—are significant determinants of their mechanical integrity, their impact has not been explicitly discussed when reporting bulk mechanical properties. Controlled architectures were designed by systematically varying porosity (30–75%, H/D ratio (0.5–2.0) and pore size (0.25–1.0 mm) and fabricated using fused filament fabrication technique. The influence of the three parameters on compressive mechanical properties—apparent elastic modulus Eapp, bulk yield stress σy and yield strain εy—were investigated through a multiple linear regression analysis. H/D ratio and porosity exhibited strong influence on the mechanical behavior, resulting in variations in mean Eapp of 60% and 95%, respectively. σy was comparatively less sensitive to H/D ratio over the range investigated in this study, with 15% variation in mean values. In contrast, porosity resulted in almost 100% variation in mean σy values. Pore size was not a significant factor for mechanical behavior, although it is a critical factor in the biological behavior of the scaffolds. Quantifying the influence of porosity, H/D ratio and pore size on bench-top tested bulk mechanical properties can help optimize the development of bone scaffolds from a biomechanical perspective.


Introduction
Porous scaffolds to guide and stimulate tissue growth are increasingly considered a viable option in bone tissue engineering applications. Optimal osteogenic signal expression and subsequent differentiation of cells seeded on the scaffold are influenced by physical scaffold parameters such as mean porosity, pore size and pore interconnectivity and mechanical parameters such as strength and elastic modulus of the fabricated bulk structure [1][2][3][4]. Porosity and interconnectivity ensure migration, attachment proliferation and differentiation of cells in the scaffold and flow for nutrient transport and waste evacuation [5]. Similarly, scaffold macro-pore size is an important variable affecting the ability of bone scaffolds to accommodate cell ingrowth and new bone formation [5][6][7][8]. Although an ideal scaffold pore size for efficient bone regeneration has yet to be determined, studies have reported viable pore sizes ranging from 100 µm up to 1200 µm [7,[9][10][11][12][13][14][15].
Different types of scaffold architectures have been implemented over the last several decades, which can be classified according to their macro-porous configuration: single random porous domain, single regular porous domain and multi-domain porous [16]. The main limitation of single random porous domain scaffolds, for example sponge-type scaffolds, is that seeded cells cannot migrate into the interior regions of the scaffold. Additive manufacturing (AM), also known as rapid prototyping, has emerged as a powerful technique to address the limitation of single random porous domain scaffolds by creating scaffolds with a single macropore domain of regular morphology, such as orthogonal arrays of channels [17]. Hence, the lack of inter-connectivity presented in random pore structures is removed and flow of nutrients through the internal architecture can be facilitated. Another advantage of AM, in contrast to conventional and subtractive fabrication, is the fabrication of tissues, organs and medical devices with complex shapes and multiple materials [18]. Fused filament fabrication (FFF), which is based on heating thermoplastic filaments to their fusion point in order to fabricate a structure in a layer-by-layer process, is a popular AM method [19]. The resolution of FFF theoretically supports a minimum feature size of 100 µm [20]. In addition, FFF is generally inexpensive and therefore, together with the described advantages, the most commonly used polymer-based three-dimensional (3D) printing method for bone tissue scaffolds.
From a mechanical perspective, the bone scaffold structure should have sufficient mechanical strength to withstand normal physiological loading during the bone regeneration phase [21][22][23]. Furthermore, the stiffness of the scaffold must be tuned according to the mechanical properties of the surrounding tissue-i.e., to enable load-sharing conditions for optimal bone growth without overloading the nascent bone. This macro-mechanical requirement is typically assessed by conducting quasi-static compression tests on fabricated test specimens to determine elastic modulus (measure of bulk stiffness) and yield stress (bulk mechanical strength). Specimens used for compression testing are fabricated with the same architectural design parameters of porosity and pore size as the bone scaffold. However, these mechanical test specimen requirements give rise to several issues.
Firstly, the strength and stiffness are often bulk values, i.e., they are based on an assessment of bulk stress computed as overall applied compressive force over bulk cross-sectional area. The bulk cross-sectional area is based on the overall specimen footprint and typically does not account for the internal porous structure of the specimen, which significantly alters the effective cross-sectional area. Consequently, the elastic modulus is an apparent elastic modulus and can vary depending on the designed porosity or pore size. For 3D printed scaffolds, there is additional variability across specimens fabricated to achieve the same designed porosity and pore size due to the limitations of precision of the 3D printing process.
Secondly, the accuracy of compression testing results for trabecular bone and biomimetic cellular solid structures is strongly affected by the presence of end-artifacts [24]. Endartifacts can broadly be classified into two categories: specimen-platen interface conditions and structural end-artifacts [25]. End-artifacts distort results more strongly in shorter specimens compared to taller specimens. To standardize the mechanical characterization of porous scaffolds, international standards of traditional polymer based-materials have been widely adopted by several research groups [26][27][28]. For instance, the American Society for the Testing of Materials' (ASTM) ASTM D695 standard for compressive properties of rigid plastics defines the standard test specimen for strength measurements to be in the form of a prism or cylinder whose aspect-ratio, defined as height/diameter (H/D), is a minimum of 2/1 [29]. Nevertheless, scaffolds studies often report compression test results with lower H/D ratios-as low as 0.15 [30][31][32][33][34]. Although this H/D ratio may be sufficient to meet minimum requirements for the continuum assumption and is adequate for biological experiments to assess cell toxicity, proliferation and adhesion, it may result in an inaccurate characterization of mechanical property.
While studies have recognized these issues on a qualitative basis, their impact, especially from a biomechanical perspective, has neither been thoroughly quantified nor explicitly discussed when reporting mechanical properties of bone scaffolds. Hence, a better quantitative understanding and awareness of the influence of porosity, pore size and H/D-ratio on bone scaffold mechanical properties is needed to optimize the development of these scaffolds for tissue engineering. Therefore, in this study, controlled bone scaffold architectures were designed by systematically varying three parameters: porosity, H/D ratio and pore size, and 3D printed with FFF. The influence of the three parameters and their interactions on scaffold mechanical properties such as Elastic modulus, Yield stress and Yield strain were investigated through a multiple linear regression adjustment by a stepwise multiple linear regression model.

Materials and Methods
Porous mechanical test specimens were fabricated with a commercial 3D printer in deliberate combinations of pore size, porosity and H/D ratio. To assess the impact of the parameters on the mechanical properties-apparent elastic modulus, yield stress and yield strain-a stepwise multiple linear regression model-based study was conducted.

Material and 3D Printing of Scaffold Test Specimens
Mechanical test specimens were fabricated with commercially acquired polylactic acid (PLA) 1.75 mm diameter filament using a desktop FFF 3D printer (Mbot Grid II+, Hangzhou, China) at 210 • C and 60 mm/s printing speed [35]. The printing parameters are listed in Table 1.

Scaffold Test Specimen Design and Fabrication Process
Mechanical test geometries were based on common scaffold designs to compare with published studies and followed the ASTM D695 standard for mechanical characterization of polymers. The mechanical testing specimens were cylindrical with a constant diameter, D, of 10 mm. Inner architectures were designed following the procedure delineated in Figure 1 by varying three main parameters: (1) Height (H/D ratio), (2) Porosity and (3) Pore size. Each parameter had three levels: low, medium and high, as explained below:

1.
Height, H: A "low" height value of 5 mm represented a 1/2 H/D ratio (D = 10 mm) and is commonly used in biological assessment. End-effects, as defined by St. Venant's principle, tend to be significant in these geometries. To minimize influence of endeffects, ASTM D695 defines an H/D ratio of 2/1. Accordingly, a "high" height value of 20 mm was defined in this study. To effectively compare the mechanical behavior and the influence of the end effects, a "medium" height value of 10 mm was additionally defined representing an H/D ratio of 1/1. Thus, the respective H/D ratios were 0.5, 1.0 and 2.0.

2.
Porosity: Scaffold designs generally mimic the porosity of bone tissue. Low porosity structures such as cortical bone range between 5-30% porosity, while cancellous bone porosity is mostly in the range of 75-95% [36]. In this study, a "low" porosity level close to 30% and a "high" level near 75% were defined. The "medium" porosity was 50%.

3.
Pore size: In the current study, pore sizes from 0.25 to 0.5 mm were defined as the "low" level. Pores from 0.5 mm until 0.75 mm were "medium" level and pore sizes from 0.75 to 1.00 mm were the "high" level. Figure 2 summarizes the different combinations of pore size and porosity for a representative specimen height of 10 mm. The designed scaffolds were printed for each condition in the horizontal printing plane, where orientation of the fibers and their bonding was enhanced over other planes [37]. Table 2 summarizes the experimental design with the scaffold design combinations. Based on the combination of parameters-six specimens each with three levels for each of the three factors (6 × 3 3 )-a total of 162 specimens were fabricated. Figure 1. Scaffold test specimen geometrical design process to match theoretical values of porosity, heights and pore size: (A) the first layer was designed, the region of interest (ROI) of the struts that are created with a width (SW), length (SH) and a height (PH), as denoted by red lines. PS corresponds to the pore size and is equal to the pore height (PH). (B) A second layer is added by rotating the first one by 90 • and placing it on top of it, a circumference with diameter (CD) is designed and everything outside it is removed producing (C). (D-F) The remained part is duplicated along the cylindrical principal axis (z-axis) as required for the specimen height. The final specimen geometry with length of 5, 10 or 20 mm was exported for 3D printing as a STL file.

Morphology Characterization
Diameter and height of each printed specimen were recorded as a mean of three measurements for each dimension, as measured with a set of calipers. Scaffold porosity was measured with a buoyancy scale following the Archimedes method [38]. In this case ethanol, with density of 0.789 g/mL was used as the liquid with known density at room temperature of 22 • C. Porosity was calculated based on the formula: where: W d is the dry weight measured before the immersion; W s the submerged weight acquired in the balance; V the overall volume; and ρ the porosity of the displaced liquid. The porosity, as calculated based on Equation (1), was the experimentally measured porosity of the 3D-printed specimens. The theoretical, design porosity-i.e., either 30%, 50% or 75%-was confirmed based on the CAD model as the effective volume of the scaffold material (total volume of the struts) divided by the bulk volume (H × π × D/4) of the cylinder. Differences between the experimentally measured porosity and the design porosity were then expressed as percentages.
Following porosity measurements, specimens were dried and stored for subsequent evaluations. Optical measurements were performed to measure the specimen pore size with digital pictures acquired by an optical microscope (Leica, Leica Camera AG, Wetzlar, Germany). Furthermore, micro computed tomography (micro-CT) scans were conducted to verify the inner structure of the samples (Figure 3). Images were acquired in an EasyTom micro (Rx Solutions, Boynton Beach, FL, USA) using a configuration of voltage of 90 kV and current of 200 µA, frame rate of 2 fps. Each scan of 360 • took 20 min to achieve a resolution of 10 µm. First, the software X-Act (Rx Solutions, Boynton Beach, FL, USA) was used to preprocess the images to generate a dataset of layers along the z axis of the scan volume. These images were analyzed in VG Studio and compared with the designed CAD for printed irregularities. The samples were imaged by placing them in a low-density material to avoid undesired rotation of the specimen while scanning.

Scaffold Test Specimen Mechanical Property Characterization
Compression tests were performed at room temperature on a universal materials testing machine (Test Resources, Shakopee, MN, USA) at a fixed, quasi-static speed of 1.27 mm/min following the standard ASTM D695 [39]. Specimens were placed in the center of the plate and a preload of 50 N was applied. The test was conducted until specimen nominal strain was at least 30% strain. Bulk stress and strain were computed as: where: F = force applied at the crosshead; CSA bulk = nominal cross-sectional area; δ = crosshead displacement; L o = initial length (height) of the specimen. Hooke's Law of elasticity in elastic solids was applied to calculate apparent elastic modulus as follows: σ = E app ×ε; where σ is the bulk compressive stress, E app the apparent elastic modulus and ε the bulk strain. Apparent elastic modulus, E app , was found by linear regression of stress-strain data from the linear segment of the test data, gener-ally between 0 and 2% strain. The Yield stress and Yield strain were determined with a 1% offset strain [40].

Data Analysis
Statistical analysis was performed using R (version 3.6.3) [41]. Results were expressed as means and standard deviations and, in all cases, the level of significance was set at α = 0.05. First, Spearman's coefficients [42] were calculated to determine the correlation between the response variables and the possible explanatory variables (Height, Porosity and Pore Size). Next, effects of Height, Porosity and Pore Size were assessed based on a Mann-Whitney-Wilcoxon Test to identify statistically significant effects on the response variables. Based on the results obtained in the Mann-Whitney-Wilcoxon Test for each explanatory variable, a multiple linear regression model was developed by using the measured experimental values of the explanatory variables and their respective interactions. The linear model was: Y = (Height) × β 1 + (Porosity) × β 2 + (Pore size) × β 3 + (Height:Porosity) × β 4 + (Porosity:Pore size) × β 5 + Intercept (3) where Y are the response variables, namely, E app , Yield Stress and Yield Strain. β i are the coefficients of the regression model associated with the variable i, namely, Height, Porosity, Pore size and their respective interactions.
Subsequently, in order to determine which variables contributed to the multiple linear regression model, a step-wise regression algorithm by the forward method [43] was applied to define the influence of the independent variables. Normality assumptions inherent to the multiple linear regression model were verified with the Kolmogorov-Smirnov test [44].

Scaffold Test Specimen Morphology Characterization
The 3D printed scaffold test specimens had consistent and uniform bulk dimensions (height and diameter), with low standard deviations and errors. The variation in height (H) across the three H/D ratio groups was ≤2%, while the variation in diameter (D) was less than 6%. Given the low variation in the bulk dimensions, H/D ratio was maintained as a categorical variable with three levels (H/D = 0.5, 1.0 and 2.0) for the statistical analysis. The variability between design and printed structures was verified by superimposing the micro CT-based volume of the printed scaffold onto the CAD model. (Figure 3).
In contrast, variation in the internal architecture of the printed samples, represented by porosity and pore size, was much larger than was found for the bulk dimensions. For example, compared to low and medium porosity, specimens with high porosity (75%) had a larger error, with 10-50% percent errors between the measured and theoretical porosity. The medium porosity (50%) samples had smaller errors, 11-21%. Samples with low porosity (30%) generally had the lowest error, 0-8%, except for an atypical error of 27% for the medium pore size sub-group (0.75 mm). The "low" pore size sub-group (0.50 mm) had the biggest percent error of 46-50%, followed by the "medium" size (0.75 mm) with 19-20% and the "large" size (1 mm), 10-11%. As a result, porosity and pore size were treated as continuous variables in subsequent statistical analyses.

Scaffold Test Specimen Mechanical Property Characterization
Apparent elastic modulus, E app , was positively correlated with specimen H/D ratio and 20 mm height specimens (largest H/D ratio = 2.0) had, on average, the largest E app , which decreased progressively for the 10 mm and 5 mm height groups. Yield strain, on the other hand, exhibited a strong negative correlation. Specimen H/D ratio did not influence the bulk Yield stress (Table 3, Figure 4D).
Pore size, as an independent variable, did not have a significant effect on E app of the specimens in this configuration; however, porosity did have an effect. At the highest levels of porosity, E app decreased, as expected with porous structures. Yield stress values also decreased with increase in porosity, while E app and yield stress were strongly negatively correlated with porosity. Yield strain exhibited a mild negative correlation ( Figure 4C).
Apparent elastic modulus and yield stress were related to porosity with an exponential decay, ae bx ( Figure 4B,D, respectively). Yield strain was linearly related to porosity ( Figure 4C). The normalized modulus (apparent elastic modulus divided by the material elastic modulus versus porosity curves) are overlayed with published curves in Figure 5 according to ASTM standard D696. A summary of the measurements from the morphology and mechanical property characterization can be found in Table 3.

Multiple Linear Regression Analysis
Statistical analyses were performed to evaluate the correlation between the different parameters. Spearman correlation was used to obtain the nonparametric measure of rank correlation. This correlation describes how well the relationship can be defined using a monotonic function. Table 4 shows Spearman correlation coefficients between the independent variables (Height, Porosity and Pore Size) and the response variables (E app , yield stress and yield strain). Out of the three response variables, height exhibited the highest correlation with yield strain (−0.761) followed by E app (0.501) and negligible correlation with yield stress (0.043). Porosity was highly correlated with both E app (−0.859) and yield Stress (−0.912), but had a low correlation with yield strain (−0.269). Finally, pore size was not strongly correlated with any of the three response variables, the highest coefficient being −0.204 for yield stress. Actual pore size was considered as a continuous variable in the statistical analysis. Apparent elastic modulus (E app ), yield stress (σ) and yield strain (ε y ) were determined from compression test data. Measurements from 3D printed scaffold test specimens were based on six replicates and results are expressed as mean ± standard deviation.     Table 5 shows the step-wise statistical analysis results, with models for the response variables (E app , yield stress and yield strain) based on the specimen parameters (height, porosity and pore size). The model successfully explained up to 96% of the variation in both apparent elastic modulus (E app ) and yield stress. For E app, porosity was the principal parameter (R 2 = 73%), followed by height (R 2 = 19%). For yield stress, the principal parameter was porosity (R 2 = 94%), which explained almost all the variation in yield stress. Finally, the model explained only up to 72% of the variation in yield strain, which was mainly represented by the height (R 2 = 60%), with porosity accounting for the remaining 12%. Table 5.
Step-wise statistical analysis with models for the three response variables. Coefficient of determination (R 2 ), coefficient of the regression model associated with the variable (β i ) and p-value (p-value) are listed for each of the variables and their interactions, were p-value < 0.0001: ****; p-value < 0.001: ***; p-value < 0.01: **; p-value < 0.05: *; p-value >= 0.05: NA (does not contribute to the model).

Discussion
One hundred and sixty-two 3D-printed scaffold test specimens with controlled geometries were fabricated to systematically evaluate the variation in the mechanical response of 3D structures obtained based on variations in H/D ratio, porosity and macro-pore size. Combined, these parameters resulted in almost a six-fold variation in the full range of apparent elastic modulus and bulk yield stress values-from 189 MPa to 1220 MPa and from 7 MPa to 41 MPa, respectively.
Results from the statistical analysis can help us understand how the parameters tested in this study affect mechanical properties.

Elastic Modulus
In Table 5, the E app is well represented with the proposed model (R 2 of 96%) with porosity as the principal influencing parameter (R 2 of 73%), similar to findings in literature [18,[53][54][55][56][57]. The negative β i suggests that the increase of porosity reduced the stiffness of the samples with a β i of −8.92 per percentage increase in porosity.
H/D ratio had a relatively smaller influence on the E app (R 2 of 19%) within the model, but a higher sensitivity on the samples with a value of β i of 43.46. Moreover, the model also revealed an interaction or mild confounding effect between the height and porosity. For example, specimens with high porosity (negative influence on E app ), but high H/D ratio (positive effect on E app ) exhibited elastic modulus values close to specimens with low porosity and low H/D ratio ( Figure 4B). A sensitivity study was carried out on the intercept value of the model to extrapolate the response of a solid sample. The model, driven by height with null porosity and thus null pore size, showed an E app of 1034, 1251 and 1684 [MPa]-a variation of almost 60%-purely due to changes in specimen height between 5 mm (H/D = 0.5), 10 mm (H/D = 1.0) and 20 mm (H/D = 2.0), respectively. These values are consistent with those reported previously with the same configuration [20], which serves as a further validation of the model. The change in E app for each different height essentially represents the effect of the H/D ratio in the mechanical response [35,58]. Thus, although the influence of the specimen H/D ratio on the elastic modulus relative to specimen porosity may be smaller, it is still significant and must be taken into account while comparing across studies.
Finally, neither the range of pore sizes nor its interaction with the other variables (height and porosity) had a significant effect on E app (p > 0.05), which is consistent with literature [18,[54][55][56][57]59].

Yield Strain
Compared to elastic modulus the relative influence of H/D ratio and porosity on yield strain were more or less reversed. The regression model had a lower predictive strength, predicting variation in the yield strain with an R 2 of 72%. H/D ratio was the principal influencing parameter (R 2 of 60%). The influence of porosity was relatively smaller (R 2 of 12%). Although the influence of these parameters on yield strain has been noted in past studies [35,58], this study quantifies these effects in a model where the sensitivity is minimum for both parameters, with β i = -0.0014 for height and β i = -0.000276 for porosity. Hence, for a given material, neither the inner architecture, nor the H/D ratio changes affected the yield strain substantially, with a consistent value of 0.07. Furthermore, the pore size did not have a significant influence on the mechanical responses of the proposed model, as has also been shown in the literature [3].

Yield Stress
Finally, similar to apparent elastic modulus, yield stress and its variations are also represented well by the model (R 2 = 96%). However, porosity variations could explain most of the variations in yield stress (R 2 of 94%) and with a low representation by the height (R 2 of 2%). The β i shows that yield stress was inversely correlated with these parameters, a trend consistent with previous studies [18,[55][56][57]59]. Notably, however, yield stress was the only response variable, where pore size displayed a significant influence (p < 0.05); nevertheless, the influence within the model was low (R 2 < 1%).
The yield stress findings can be understood as essentially a product of E app and yield strain. Further, E app is strongly correlated with yield stress, while it is poorly correlated with yield strain (Figure 6). Both E app and yield strain are strongly negatively influenced by porosity, resulting in an extremely strong negative effect of porosity on yield stress. On the other hand, while E app was positively affected by specimen H/D ratio ( Figure 7A), yield strain was negatively affected ( Figure 7B), essentially cancelling out the effect of specimen height on the yield stress ( Figure 7C). Consequently, yield stress appeared less sensitive to variations in specimen H/D ratio.  . IQR = interquartile range between first to third quartile and "n" = number of non-missing observations within the group. Non-overlapping notches represent significant differences [60,61].

Conclusions
Within the context of the three factors investigated together in the current study, H/D ratio and porosity of the fabricated structures had a strong influence on the mechanical properties commonly used to understand the mechanical behavior of these structures, namely apparent elastic modulus, yield strain and yield stress. Thus, when comparing across studies, for example, between printing techniques or even choosing candidate polymer materials for printing scaffolds, it is important to note the differences in the H/D ratios as well as porosities of the specimens used in the respective studies. Particularly for porosity, the variations in actual porosity of the fabricated specimens with respect to the designed value may also be significant enough to influence the reported mechanical properties. Depending on the specific mechanical property in consideration, either porosity or H/D ratio may have a dominating influence on the results; nevertheless, variations in both should be taken into account. Although H/D ratio appeared to significantly influence the stiffness (elastic modulus) and yield strain, yield stress did not seem sensitive to this factor within the specific range of H/D ratios investigated in this study. Thus, yield stress could potentially be a benchmark mechanical property for comparisons across studies for specimens with different heights, as the height does not have a high effect on yield stress. On the other hand, if the samples have different porosity but the same height, the yield strain is a suitable result variable for comparisons. Funding: This research was funded by "Agencia Nacional de Investigación y Desarrollo" (ANID) from the Government of Chile, through the Research Project FONDECYT Inicio" No. 11170957. J.I.C.R. thanks "Universidad Adolfo Ibáñez" for providing funding for doctoral studies by "Beca de Excelencia FIC-EMPA".
Institutional Review Board Statement: Not applicable.

Data Availability Statement:
The data presented in this study can be made available on request from the corresponding author.