Fluid dynamic parameters of naturally derived hydroxyapatite scaffolds for in vitro studies of bone cells

Hydroxyapatite scaffolds obtained from the biomorphic transformation of wood are characterized from a fluid dynamic point of view. Such material of recent introduction offers great advantages for the in vitro study of bone cells, mostly in virtue of its peculiar porous structure. Determining the flow resistance and morphological parameters of these scaffolds is an essential step towards their practical use in bioreactors and microfluidic devices. To this aim, a series of tests involving a draining fluid are performed on a set of disc-shaped scaffolds, followed by the microscopy analysis of the pores visible on the sample faces. Contrarily to what expected, a temperature dependence is observed for the flow resistance, even after normalizing it by the fluid properties. The interpretation of the experimental results is assisted by numerical outcomes from Computational Fluid Dynamics modelling, which underline some limitations in the application of classical laws to the present problem. While the complex and variable internal structure of the scaffolds prevents the systematic use of simplified formulae, a correlation is found between the flow resistance and the pore geometry, which can facilitate the characterization of further samples.


Introduction
Understanding the mechanisms governing bone transformation and healing is at the basis of bone regenerative medicine.To this aim, it is of great importance to try to re-create in vitro the same environment found inside living organisms, including the flow conditions.In fact, it is well ascertained that bone cells, i.e. osteoblasts, osteoclasts and osteocytes, regulate the process of bone growth and resorption in response to external mechanical stimuli, which are sensed in the form of fluid shear stresses (FSS) acting on their surface [1].These are in the order of the Pascal or lower, depending on the cell type and on the physiological conditions [2].Fluidic systems thus need to be designed to be able to generate the required levels of FSS.
One of the key elements to be integrated in culturing systems is the so-called scaffold, that is, an extracellular matrix providing the attachment sites for the cells.Three-dimensional scaffolds are considered superior to planar scaffolds, because they allow for cell expansion and re-arrangement in all directions, but their fluid dynamic characterization is not always straightforward.In the current investigation, hydroxyapatite scaffolds resulting from the biomorphic transformation of wood are adopted as the substrate for bone cells to be cultured under perfusion.This naturally derived material, termed as B-HA, has been recently developed [3] and presents a 3D structure with interconnected, hierarchically organized pores that closely resembles the microstructure of real bone.Its inherent flow resistance is an essential parameter that must be determined in order to calculate and control the flow rate and the shear stresses the cells inside the pores are exposed to.However, as the result of the specific fabrication process, the internal geometry of B-HA is anisotropic, not predictable and specific of the individual scaffolds.Furthermore, it is hard to be directly inspected, since B-HA is prone to shatter upon sectioning with common laboratory instrumentation.
Since those factors prevent the application of simple formulae, a series of experiments is carried out here to accurately derive the fluid dynamic parameters of a set of B-HA samples.The flow resistance is reconstructed by exposing the scaffolds to a draining volume of fluid and the scaffold morphology is analysed via microscope observations.The collected data on pore number and size are used in combination with the values relative to flow resistance to extract useful correlations in order to simplify the characterization of B-HA.The experimental campaign is supported by a Computational Fluid Dynamics (CFD) study, aimed at investigating in detail the flow development inside the pores.The numerical outcomes obtained on one of the available samples are compared with the experimental results and used to explain the observed trends.

Experimental methods
The naturally-derived B-HA scaffolds adopted in the present study are in the form of round pills nearly 10-11 mm in diameter and 4-5 mm in height.To obtain accurate values of the flow resistance, perfusion tests are performed on a set of 12 scaffolds.As depicted in Figure 1a, each sample is wedged in the orifice of a hourglass-type transparent container and exposed to a draining volume of bi-distilled water (nearly 500 ml).An imaging system is used to capture the temporal variation of the liquid level in time lapse mode.The acquired frames are analysed by means of an image processing software (ImageJ, NIH, USA) to extract the instantaneous liquid height.More details on the experimental apparatus and data elaboration can be found in [4].The scaffold flow resistance, Rf, is then computed for each frame from the following formula, after integration between the starting time and the considered instant: In the above equation, Δp is the pressure difference over the scaffold dictated by the hydrostatic head, being ρ the fluid density and h the hydraulic head with respect to the scaffold.The volumetric flow rate, Qv, is expressed in terms of the descending liquid level and the diameter of the liquid free surface, D, which is a function of h too (the portion of interest of the container has a conical geometry).The symbols used, together with the relative definition and units, are reported in Table 1.An average value of Rf is then extracted for the whole sequence with the associated standard deviation as the uncertainty.The experiment is repeated 4-5 times both at room temperature and inside an oven set at 37 °C, i.e. the temperature typically imposed during in vitro studies of cell cultures.In this way, any possible temperature dependence of the scaffold flow resistance is taken into account.The actual temperature of the liquid is measured with a thermocouple immediately before and at the end of each test.Information on the pore number and size is collected by inspecting the top and bottom surfaces of the samples with a stereomicroscope (Nikon SMZ800, Plan 1X objective).It is observed that fluid immersion induces swelling in the B-HA solid matrix, causing a general contraction of the pores which affects the scaffold flow resistance.Hence, the morphological analysis of each sample is postponed until the completion of the relative permeation tests in order to get more representative geometrical data.Microscope images are processed with ImageJ: image segmentation is performed first via filtering, rolling-ball background subtraction and semi-automatic thresholding to separate the features of interest (i.e. the pores) from the background.The obtained binary image is then scanned with the automated particle detection algorithm included in the software, which counts and measures the cross-sectional area of the detected pores.Due to the presence of pores significantly different in size, the overall procedure is repeated twice or three-times on each image by adjusting threshold as well as the particle size range in the detection algorithm.In this way, both large and tiny pores, as small as 47 μm 2 in area, are detected while reducing area over-/under-estimation and duplicated counts.In Figure 1b, the output of the processing technique applied to one of the images is displayed as an example.Eventually, pore cross-sections are approximated as circles, and the equivalent diameters are computed.Results are grouped in 11 classes of increasing width to have a proper representation of diameter distribution.

Experimental results
The thermo-physical properties of culture media commonly adopted for in vitro cell growth are slightly different from those of water [5].As a consequence, the measured scaffold flow resistances should be conveniently made independent of the test liquid properties.By supposing laminar flow conditions (as confirmed below), the Darcy's law can be considered valid on first approximation for the present problem, reading as: From a comparison with equation ( 1), it can be seen that the fraction on the right hand of equation ( 2) must be equal to Rf, as given by the dynamic viscosity of the liquid (μ) and by the height (L), crosssectional area (A) and permeability of the scaffold (k).Then, flow resistance information related solely to scaffold properties can be obtained if scaling Rf by the water viscosity at the corresponding test temperature.The Rf/μ so derived for the investigated samples distribute over one order of magnitude, ranging from (1.51 ± 0.06) •10 11 to (1.254 ± 0.007) •10 12 m -3 , relatively to room temperature conditions (21-23 °C).Such a wide interval of data is symptomatic of a great variability of B-HA internal structure, directly connected to its natural origin.
Contrarily to what expected, the Rf/μ values do not remain fixed when moving to higher temperatures (35-38 °C), though the variation of water viscosity could be considered the main responsible for any temperature dependence.A general increase is observed, with the shifted range going from (1.83 ± 0.02) •10 11 to (1.53 ± 0.04) •10 12 m -3 .This cannot be attributed to a swelling of the scaffold due to hygroscopic effects, since the two test conditions are not always applied in the same order.In Figure 2a, the relative variation of Rf/μ with respect to the room temperature value is reported for each sample, normalized by the specific temperature difference.A total mean increment of (1.08 ± 0.13) •10 -2 °C-1 is estimated when passing from an average liquid temperature of 22 °C to 37 °C.Data collected from stereomicroscope images confirm that B-HA samples present a wide variety of pore distributions, as exemplified by the bar charts shown in Figure 2b, which are relative to the two samples generating the extreme values of Rf/μ.However, as described in [6], also samples giving comparable flow resistances can exhibit a significantly dissimilar morphology.In general, the number of pores detected on the scaffold faces varies from nearly 400 to more than 900, with minimum diameters of 8 μm observable in all the cases.The maximum pore diameter, instead, differs from sample to sample, ranging from nearly 260 to 910 μm.Interestingly, the pore distribution on one face of the scaffolds roughly replicates the opposite one, suggesting that pores extend preferentially along the vertical direction.This allows for the samples to be approximated to some extent as a network of cylindrical channels.

Numerical model
In order to support the interpretation of the experimental results and to inspect in detail the flow development inside the pores, numerical simulations of the perfusion tests are performed.A 3D model is created of one of the samples, having intermediate fluid dynamic and morphological characteristics within the set of the measured scaffolds.The pore diameter distributions extracted from the front and rear face of the scaffold are averaged, and the resulting distribution (Figure 3a) is used to populate the model with cylindrical pores of length equal to the scaffold height (4.3 mm).The measured scaffold diameter (10.40 mm) is preserved too.The model is built by imposing circumferential symmetry, so that a 30° slice, like the one sketched in Figure 3b, hosts the same relative pore distribution than the full sample.A 1 mm thick fluid plenum is added above and below the scaffold to provide numerical stability and to reproduce the spontaneous flow subdivision into channels of different size, according to the individual channel resistance.The computational tool SIMCENTER STAR-CCM+ (SIEMENS SISW, v. 2021.3),based on the Finite Volume method, is used to solve the continuity and Navier-Stokes equations numerically.The experimental configurations are replicated by specifying thermo-physical properties of water at two different temperatures, namely 22 and 37 °C, which correspond to the average liquid temperature in flow resistance tests at room temperature and in an oven, respectively, on the selected sample.For simplicity, the liquid level is kept fixed, associated to a mean pressure difference of Δp = 947 Pa over the scaffold.This implies a Reynolds number inside the pores not greater than Re = 80, as computed from the mean experimental data.Hence, laminar flow conditions are applied while adopting a steadystate solver and second-order schemes for discretization.
By taking advantage of the model axial symmetry and the preferential direction of the flow, only a 30° slice is meshed and solved, allowing for a reduction of the computational requirements.An unstructured quadrilateral surface mesh is generated separately for the scaffold/channel region and the two plenums, with a minimum cell thickness of 1 μm at the walls of the smallest pores.Both grid types are axially extruded with a uniform spacing of 0.1 and 0.05 mm, respectively, for a total of 1.1 M cells.In Figure 3c, a detail of the mesh in the channel region can be observed.Stagnation and static pressure are imposed at the inlet of the upper plenum and at the outlet of the lower plenum, respectively, whereas a no-slip condition is specified at the liquid/solid boundaries.The reader is referred to [6] for further details on the numerical model, including the assessment of the mesh quality.

Numerical results
Normalized flow resistances are extracted from the simulation outputs, by dividing the imposed pressure difference by the obtained volumetric flow rate and by the temperature-dependent water viscosity.Numerical and experimental results for the selected sample are compared in Figure 4a.It can be noted that CFD predictions tend to overestimate Rf/μ by 20-30%, possibly because the boundary condition specified at the liquid/solid interface does not adequately account for the absorption mechanisms that are responsible for the swelling of the scaffolds.Dedicated investigations are needed in this sense for future improvements of the numerical model.Nevertheless, the increase of Rf/μ with temperature is still evident.
To better elucidate this point, contours of the velocity magnitude plotted at the two temperatures are depicted in Figures 4b and 4c.Higher flow velocities are associated with the highest temperature, following a reduction in water viscosity.More noticeably, a non-negligible entry region, extending over 1/4 of the channel length, can be seen at the channels inlet already in the case at the lowest temperature, and the temperature increment induces its further expansion.The absence of a fully-developed flow all along the pores poses some limitations to the applicability of the Darcy's law to the present problem, and justifies the preservation of the temperature dependence in the Rf/μ values.

Discussion
Under the approximation of parallel cylindrical pores having the same length, the geometrical data acquired from the microscope images could be used to estimate an equivalent normalized flow resistance (Rf eq/μ), according to the following formula: Rf i is the resistance of the i-th pore having diameter Di, and is calculated from the Hagen-Poiseuille's equation.For each sample, equation ( 3) is applied separately to the pores detected on the front and rear face, and results are averaged.However, the vast majority of resistances obtained in this way seem to underestimate the actual values, showing discrepancies up to 87% relatively to the measured ones, as demonstrated by the data appended in Table 2, where measured data referred to tests in oven are reported.One of the causes of such differences is identified in the limited field of view of the stereomicroscope in use, which does not allow for the detection of the most peripheral pores even at the lowest magnification.Though pores near the specimen borders are less abundant, as shown in [4], and do not significantly affect average outputs, their contribution in determining the equivalent flow resistance is lost.To fix this point, a microscope mounting a lens with a slightly lower magnification would be needed.On the other hand, less resolved images would impair the detectability of the smallest pores, and would aggravate some errors already introduced by the image processing technique.These include, for example, the interpretation of adjacent pores as a single, merged structure, which implies a reduction of the value calculated with equation (3).Another explanation for the observed discrepancies comes from the complex internal structure of B-HA, which in some cases can invalidate the geometrical simplifications made.
In light of the above considerations, perfusion tests stand out as the best method for an accurate determination of B-HA flow resistance.However, this appears as an unfeasible solution when multiple scaffolds need to be analysed.Introducing appropriate correlations between the fluid dynamic and statistical morphological properties of the scaffold would represent a helpful alternative.For this purpose, flow resistances are further normalized in order to get non-dimensional quantities.In particular, Darcy number (DaL) is computed referred to the scaffold height, as given by equation ( 4):  Since the operating temperature in typical scaffold applications is T = 37 °C and temperature dependence has been observed for Rf/μ, resistance values in equation ( 4) are inserted pertaining to tests carried on in oven only.
On the other hand, different non-dimensional statistical parameters describing pore geometry are considered, including: (i) the mean pore diameter scaled by the scaffold diameter, (ii) the mean pore area (Amean) scaled by the captured scaffold area (Asc), and (iii) surface porosity (NAmean/Asc, being N the number of pores).For their determination, averaged data from the front and rear face are used for each sample.It is pointed out that the scaffold area entering the surface porosity and the Amean/Asc ratio is only the visible portion of the scaffold cross-sectional area in the microscope images.Moreover, surface porosity and scaled pore diameter/area are parameters related to the scaffold surface.Their involvement as the independent variable in correlations is reasonable only for samples which have an approximately regular internal structure.For this reason, samples exhibiting Rf eq/μ largely different from the measured values (> 70%) are excluded in the search of correlations.Among the explored possibilities, the best correlation is found for DaL as a quadratic function of Amean/Asc, which is reported in Figure 5.The achieved relation returns flow resistance values that, in the vast majority of cases, are more reliable than those resulting from equation ( 3), but still have some discrepancies with respect to the measured ones.Hence, while using this function would allow for a sensible reduction of experimental and analysis efforts in the extraction of B-HA fluid dynamic properties, its validity needs to be evaluated by considering the scaffold final application.
Figure 5. Plot of the Darcy number as a function of the normalized pore mean area.The reported correlation is found from green data ("selected data") only.

Conclusions
The flow resistance and pore geometrical parameters of naturally derived hydroxyapatite scaffolds are determined as functional parameters to control the stimuli over bone cells cultured under dynamic conditions.The wide scattering of the extracted data directly reflect the variable internal structure of this type of scaffold.For modelling purposes, this can be approximated to a network of cylindrical pores, simplifying the investigation of flow development by CFD analysis.The inspection of velocity contours from numerical simulations reveals the presence a non-negligible entry region inside the pores, which undermines the validity of the classical Darcy's law when decoupling flow resistance from the fluid properties.On the other hand, for predictive purposes, the simplification of the internal pore geometry can lead to quite unreliable results.A correlation between the fluid dynamic properties of the scaffolds at their typical operating temperature and pore statistical parameters is proposed as a better alternative.The correlation found still suffers from some limitations, consequent to the assumptions made, but can be reasonably employed to streamline the preparatory characterization of the scaffolds to be used in in vitro studies.

Figure 1 .
Figure 1.a) Schematic of the test device and microscope image of one of the scaffolds (upper face) inserted into the device orifice.b) Pores identified on the scaffold face (detailed view) via repeated processing procedure.

Figure 2 .
Figure 2. a) Relative variation of scaffold flow resistance with respect to the room temperature value divided by the average temperature difference between the two test conditions.b) Pore distribution for two scaffolds having relatively low (Sample 10) and high (Sample 6) flow resistance.

Figure 3 .
Figure 3. a) Average pore distribution of the sample selected for the numerical analysis.b) Arrangement of the pores in a 30° slice of the model (top view).c) Detail of the grid generated in the scaffold region.

Figure 4 .
Figure 4. a) Comparison of the experimental and numerical results for Rf/μ for the selected sample.b) and c) Contours of the velocity magnitude for two different temperatures imposed to the water properties: b) T = 22°C and c) T = 37 °C.

Table 1 .
List of symbols.

Table 2 .
Comparison of the normalized flow resistances for the set of scaffolds: values measured from the tests in oven (Rf ov /μ) against values obtained from equation (3) (Rf eq /μ).