Analysis of Gene Expression Signatures for Osteogenic 3d Perfusion-bioreactor Cell Cultures Based on a Multifactorial Doe Approach

The use of multifactorial design of experiments (DoE) in tissue engineering bioprocess development will contribute to the robust manufacturing of tissue engineered constructs by linking their quality characteristics to bioprocess operating parameters. In this work, perfusion bioreactors were used for the in vitro culture and osteogenic differentiation of human periosteum-derived cells (hPDCs) seeded on three-dimensional titanium (Ti) alloy scaffolds. A CaP-supplemented medium was used to induce differentiation of the cultured hPDCs. A two-level, three-factor fractional factorial design 640 was employed to evaluate a range of bioreactor operating conditions by changing the levels of the following parameters: flow rate (0.5–2 mL/min), cell culture duration (7–21 days) and cell seeding density (1.5 × 10 3 –3 × 10 3 cells/cm 2). This approach allowed for evaluating the individual impact of the aforementioned process parameters upon a range of genes that are related to the osteogenic lineage, such as collagen type I, alkaline phosphatase, osterix, osteopontin and osteocalcin. Furthermore, by overlaying gene-specific response surfaces, an integrated operating process space was highlighted within which predetermined values of the six genes of interest (i.e., gene signature) could be minimally met over the course of the bioreactor culture time.


Introduction
Bone tissue-engineering (TE) strategies aim at the ex vivo dynamic cultivation of human mesenchymal stem cells, seeded on appropriate three-dimensional (3D) biomaterials (scaffolds), and the use of bioactive factors within a bioreactor setup in order to create TE constructs with regenerative potential upon in vivo implantation [1].Scaffolds allow cell attachment and proliferation within their pores, while at the same time providing the surface properties and porous structure that are required to permit cell growth in a setup that resembles more the native tissue than 2D culture conditions [2,3].In order to support improved nutrient and oxygen mass transport to cells residing at the center and not only at the peripheral layers of the scaffold, while at the same time more actively controlling the cellular microenvironment, TE constructs are often cultured in (perfusion) bioreactors [4][5][6].Furthermore, with the use of bioreactors for TE construct culture, critical culture parameters, such as dissolved oxygen concentration, pH and perfusion velocity, as well as the magnitude of the applied shear stress to the cells can be controlled [7].
However, currently, no quantitative relationships exist that link cell phenotype during bioreactor culture to operating bioprocess parameters [8].The clinical translation of in vitro produced TE constructs as graft material requires that process parameters are determined and controlled in such a way that they enable cellular constructs to be robustly manufactured [9].The selected process parameters will ultimately determine the type, quality and quantity of the produced tissue, which is typically assessed as a function of numerous output variables or performance endpoints (e.g., cell content and/or gene expression levels).In a subsequent phase, these process parameters could be used to actively steer the properties of the in vitro-developed TE construct to mimic those of the target tissue, which, in this study, is "bone" [10].To this end, the design of experiment (DoE) approaches, with a track record in bioprocess optimization, have been introduced for stem cell [11][12][13] and tissue engineering [14,15] research.In these DoE studies, multiple-factor screening of stem cell/TE process parameters and their ranking in terms of significance were obtained.Furthermore, a better understanding was gained of the possible interactions of these process parameters, which, in turn, will help in deciphering the complex relationships that are entailed when stem cell populations interact with their culture environment [16].
A number of studies have shown that during bioreactor culture of osteoprogenitor cells, shear stress resulting from the perfusion of culture medium not only improves the spatial distribution of the cells in the scaffolds, but that it also promotes the mineralization of the secreted extracellular matrix [17,18].Furthermore, shear stress induces cell differentiation towards the osteogenic lineage indicated by the increased expression of early osteogenic markers, such as, collagen type I (Col1), alkaline phosphatase (ALP), transcription factors like osterix (OSX) and late differentiation markers, such as osteopontin (OPN) and osteocalcin (OCN) [19][20][21].However, it is expected that gene expression signatures under perfusion conditions might be application specific, as they are dependent on the presence of various bio-instructive factors in the culture medium, like, for example, Ca 2+ and P i ions [22].In addition, certain culture parameters might be interdependent [16]; the cell seeding density influences the spatial distribution of cells [15], their subsequent 3D growth [23] and also their differentiation [24,25], but it also correlates to the amount of secreted growth factors in the medium and might therefore mediate the expression of osteogenic markers.This makes the interpretation of these gene expression signatures difficult and results in a limited understanding of the osteogenic differentiation process of cell-based constructs in perfusion culture.One of the key objectives in the field of bone TE will be to efficiently guide the osteogenic differentiation process of cultured mesenchymal stem cells towards a predefined gene expression signature that would ensure the clinical applicability of the produced TE constructs [26].However, to date, there is minimal knowledge regarding the determination of those minimal requirements for in vitro produced TE constructs, in order to ensure a successful in vivo outcome.Hence, the systematic production of TE constructs in well-characterized bioprocesses is a prerequisite for compliance with regulatory guidelines and subsequently successful translation to the market.
In this work, the TE constructs of interest were human periosteum-derived cells (hPDCs) seeded on regular pore Ti6Al4V (Ti) scaffolds and cultured in an osteogenic medium.The aim of this study was to gain insight into the effect of standard culture process parameters, such as seeding density and perfusion flow rate, on the in vitro osteogenic cell differentiation as a function of culture time.Gene expression analysis of the hPDCs in TE constructs after osteogenic 3D perfusion culture was performed using a multifactorial DoE approach.This approach allowed the evaluation of multiple process parameters simultaneously in a cost and time efficient way, as opposed to the one-factor-a-time approach [27].Based on the obtained data, a proof-of-concept was given that, for a dedicated TE construct process set-up, an integrated "window of operation" could be determined that resulted in a cell population with a given gene signature along the osteogenic differentiation pathway.

Ti6Al4V Scaffolds
Porous Ti6Al4V scaffolds (Ø = 6 mm; h = 6 mm) with a highly regular design were produced by selective laser melting [28].Based on micro-CT image analysis, the following scaffold characteristics were determined: 73% ± 1% porosity, 6.5 ± 0.2 cm 2 surface area, 245 ± 2 µm strut diameter and 755 ± 3 µm pore size.Scaffolds were cleaned in an ultrasonic bath consecutively in pure acetone, 70% ethanol and distilled water for 30 min in total.Subsequently, scaffolds were oxidized for 12 h in a 5 M sodium hydroxide solution at 60 °C, rinsed with distilled water and autoclaved.One day before cell seeding, scaffolds were pre-wetted with culture medium by vacuum impregnation, incubated for 2 h in an incubator at 95% humidity and dried overnight in a non-humidified incubator.

hPDC Culture
hPDCs were isolated from periost biopsies of different donors and pooled as described by Eyckmans et al. [29].The ethics committee for Human Medical Research (KU Leuven) approved the procedure, and biopsies were taken from the patient after informed consent.The cells were cultured in a monolayer in growth medium consisting of high glucose Dulbecco's modified Eagle's medium (Invitrogen) supplemented with 10% fetal bovine serum (FBS; Lonza BioWhittaker, Basel, Switzerland), 1% sodium pyruvate and 1% antibiotic-antimycotic (100 units/mL penicillin, 100 mg/mL streptomycin and 0.25 mg/mL amphotericin B; Invitrogen).At 80%-90% confluence, the cells were trypsinized with Tryple Express (Invitrogen) and replated for up to 5 passages.At passage 6, the cells were trypsinized and seeded in the pre-wetted Ti scaffolds using static drop-seeding at a density of 10 5 or 2 × 10 5 cells per scaffold.The scaffolds were cultured simultaneously in an in-house-developed perfusion bioreactor, in which culture medium was pumped by a peristaltic pump (PC-24; Ismatec SA, Vancouver, BC, Canada) from a disposable medium reservoir (50-mL Falcon tubes; Franklin Lakes, NJ, USA) containing 10 mL of culture medium, through a scaffold that was fixed in the bioreactor chamber.TE constructs were cultured for up to 21 days in a bio-instructive medium, consisting of growth medium supplemented with 6 mM Ca 2+ , 4 mM phosphate ions (P i ) and 0.05 ng/mL ascorbic acid [22].

Multi-Factorial Analysis
A fractional factorial design with three factors (flow rate, culture time and cell-seeding density) and two levels was used.The factor levels were derived from the literature and previous experimental findings (Figure 1).The 0.5 mL/min to 2 mL/min flow rate range is deemed to be within the beneficial flow rate range for osteogenic differentiation [20].hPDCs at 7 days of culture are considered to be in an early differentiation stage, while they will be expected to be fully differentiated after 21 days [30].
The low cell-seeding density of 10 5 cells per scaffold is the minimal amount of cells that [31] supports growth in the bioreactor.The high cell-seeding density of 2 × 10 5 cells per scaffold was chosen, since this is the maximum cell amount for this specific scaffold type that results in a similar cell-seeding efficiency compared to the low seeding level [15].The gene expression response on these different factor levels was determined as an average value of triplicate measurements on constructs produced in three parallel bioreactor runs per condition.As a fractional factorial design (Resolution III, main factor effect aliased with second order interactions) is used in this study, linear behavior between process parameters and gene expression response is assumed for the process space under investigation.Gene expression results were analyzed by ANOVA, and the data was fitted to a general linear model with a statistical software package (JMP Software, version 9; SAS Institute: Cary, NC, USA).The statistical significance level was set at p < 0.05.Statistical models were accepted when there was no lack of fit and no correlation in the residual plots, and the residuals were normally distributed.General linear-model parameter estimation was performed by least-squares estimation.The main effects are included in the statistical model following, (1) where: Y is the predicted response; b is the parameter estimate; X is the coded value of the factor levels; and e is the residual error.

DNA Content Analysis
The DNA content of hPDCs on scaffolds was determined using a quantitative fluorometric DNA assay (Quant-iT™ dsDNA HS kit, Invitrogen, Carlsbad, CA, USA).Scaffold samples were rinsed with phosphate-buffered saline (PBS) and lysed in 350 µL RLT lysis buffer (Qiagen, Venlo, The Netherlands) supplemented with 3.5 µL β-mercaptoethanol (Sigma, St. Louis, MO, USA).Prior to analysis, the samples underwent a freeze-thaw cycle at −80 °C.After thawing, samples were vortexed vigorously for 1 min and spun down for 1 min at 13,000 rpm.For each sample, 10 µL of supernatant was diluted in 90 µL of Milli-Q water, after which, the DNA content was quantified with a Qubit ® Fluorometer (Invitrogen), as described by Chen et al., 2012 [32].

Live/Dead Assay
A Live/Dead viability/cytotoxicity kit (Invitrogen, Carlsbad, CA, USA) was used according to the manufacturer's instructions to qualitatively evaluate cell viability and distribution on the scaffolds.Prior to capturing images by fluorescence microscopy (Zeiss Fluorescence Stereo-Microscope, Oberkochen, Germany), the samples were rinsed with PBS and incubated at 37 °C in the fluorescent dye solution (0.5 µL calcein and 2 µL ethidium homodimer in 1 mL PBS) for 20 min.

Scanning Electron Microscopy Analysis
Cell morphology on scaffolds was assessed by scanning electron microscopy (SEM) imaging coupled with an energy dispersive X-ray (EDAX) analysis (FEI XL30 FEG, Hilsboro, OR, USA) at 10 kV.The following sample preparation was used: rinsing with PSB, fixing with 2.5% glutaraldehyde (in PBS) for 1 h and postfixing in 1% osmium tetroxide for 2 h before dehydration in 50%, 75%, 95% and 100% ethanol solutions.Finally, the samples were chemically dried with hexamethyldisilazane for 3 min and gold-sputtered before SEM analysis.

Picrosirius Red Collagen Assay
3D collagen matrix production on porous scaffolds was characterized by Picrosirius Red staining (1 mg/mL Sirius Red in saturated picric acid, n = 3).The stained samples were thoroughly washed with distilled water to remove unbound dye and dried at 37 °C before qualitative analysis by stereomicroscopy.

Alizarin Red Mineralization Assay
For each bioreactor condition, mineralization of the cell cultures was assessed by Alizarin Red staining (pH 4.2) after 7, 14 and 21 days of culture.Briefly, cell cultures were rinsed with PBS, followed by fixation in 4% formalin (in PBS) for 10 min and washed with distilled water before staining with the Alizarin Red solution for 1 min.Nonspecific staining was removed by thoroughly washing in distilled water over 24 h.

Cell Proliferation
Figure 2 shows Ti6Al4V scaffolds used for cell seeding and subsequent cell culture.Cells started bridging strut junctions early on (Day 2) in the cell culture, thus initiating a 3D cell growth.Figure 3A,B shows live cells on a representative TE construct obtained after 21 days of perfusion bioreactor culture (B4 experimental condition; see Figure 3), indicating a homogeneous distribution of living cells and ECM throughout the whole scaffold.Cells seem to have filled the scaffold voids and branched out, creating a cell-ECM network in three dimensions.SEM images (Figure 2C,D) show the presence of cells in addition to CaP grains and ECM fibers in the TE construct.In static setups, 3D cell growth has been studied recently for hPDC cultures [33] and hBMSCs [34], and a two stage in vitro cell growth was described for MC3T3-E1 pre-osteoblast cells [35].In both cases, the pore geometry was seen to play an important role on cell growth kinetics.This highlights the possibility to create 3D in vitro tissue-like structures where cells may grow independently of the initial scaffold geometry and where the interaction between the cells and their own ECM and culture environment may be studied independently.Post-culture DNA content for the selected experimental conditions (as measured via DNA content analysis) is shown in Table 1.The Pareto chart of standardized effects shows that culture time and flow rate both significantly influenced hPDC proliferation (DNA content) during bioreactor culture, whereas the initial seeding density was not seen to have a significant effect on the final DNA content.For increasing culture times, a significant cell proliferation was observed, as expected for cell culture systems.Increasing flow rates, however, resulted in a slight decrease in final DNA content.In the setup used for this work, the number of cells seeded on the scaffold did not have a significant effect on the final cell number in the TE construct, as also observed elsewhere for hMSCs [37].Similar observations showing a weak negative effect of flow rate on the DNA content have been recently reported for a 15-fold increase in flow rate [38], while significant negative effects have also been observed for increasing flow rates in other bioreactor systems [17].McCoy et al. [39] have recently investigated the detachment of cells from scaffolds during perfusion culture, which was attributed to shear stress levels and to a bridged cell-scaffold attachment morphology.This phenomenon could also explain the negative effect of high flow rates observed in this work.(B)

Osteogenic Cell Differentiation: An Overview
Genes, such as Col1, ALP, OSX, OCN and OPN, have been used as markers for evaluating the specification of stem cell and osteochondro-progenitor cell types towards the osteogenic lineage [40,41].As shown in Figure 4, the use of different culture process conditions (based on the multifactorial DoE approach) resulted in significant changes in gene expression of the aforementioned genes.Overall, it may be stated that early osteogenic markers of osteoprogenitor cells, like ALP and Col1, show a significant time-dependent decrease in expression from Day 7 to Day 21, as expected [40].From Day 14 onwards, no changes were observed anymore between the different conditions for ALP.The downregulation of Col1 observed for later time points could be explained by the fact that cells become confluent within their 3D growth configuration, hence modulating the gene expression level due to less need for matrix synthesis in such a configuration.For the late osteogenic markers, such as OCN, OPN and the transcription factor, OSX, a clear time-dependent increase in gene expression was observed, as reported elsewhere for Ca 2+ and P i containing media for 2D cultures of hBMSCs [42] and during static 3D culture of hPDCs [22,30].Osteogenic differentiation was also evidenced by collagen type I secretion and mineralization (Figure 5).However, due to the DoE approach used in this work, all process parameters were systematically changed simultaneously (as shown in Figure 1), making direct comparison of the experimental conditions meaningless; thus, an integrated analysis is further required, as described in the subsequent section.

DoE Gene Expression Analysis: Flow Rate and Seeding Density Effect
In order to further interpret the influence of the selected process parameters on the gene expression levels, response surface models were created based on least-squares linear fits.The goodness-of-fit of the statistical models with the experimental data is shown in Figure 6.The strong upregulation of OSX, OPN and OCN was captured accurately, resulting in R 2 values above 0.8.The same was seen in the case of the downregulation of Col1 (R 2 = 0.79).However for the culture times evaluated in this study, both ALP and Sox9 were not clearly affected by the process conditions, resulting in R 2 values that were substantially lower (0.53 and 0.58, respectively).This demonstrates the limitation of a fractional factorial design (and least-squares linear fitting) in accurately capturing gene expression patterns that might fluctuate over the culture period of interest.However, in those cases where the u-regulation of genes was strong, the obtained fits were satisfactory.Figure 6 shows the positive effect of an increased seeding density and flow rate, which both resulted in a significant increase in the expression of the osteogenic markers, OSX (p = 0.02), OCN (p = 0.03) and OPN (p = 0.0008).Increased seeding densities may result in confluency at earlier time points, inducing osteogenic differentiation at earlier time points and, thus, enhancing the upregulation of these genes from early on in culture; hence, the positive correlation of seeding density with upregulation of these genes [25,43].However, this positive correlation has been seen to level off after a critical value has been exceeded [44].Shear stress has been linked with the upregulation of OSX, OCN and OPN during perfusion bioreactor culture [37,[45][46][47], but only a limited number of studies used human mesenchymal cells in 3D setups [19,38,41,48].Prior studies have shown that the expression of these markers increased in the presence of medium flow relative to static cultures [49], and other studies have correlated flow and/or shear stresses with an increased induction of osteogenic differentiation and tissue formation [18,50].
These findings agree with our results on the significant effect of flow on the aforementioned genes.However, we should note that, due to differences in scaffold material and design and the composition of osteogenic media, the extent and timing of osteogenic differentiation would be expected to vary across systems.Regarding ALP gene expression, our data show a negative correlation to flow rate magnitude (between seven and 21 days of culture), which has also been observed by other groups.Higher shear stress in a perfusion bioreactor decreased ALP expression after 14 days of culture [21], while a decrease in ALP specific activity from Day 7 to Day 21 was also observed during 3D culture of hMSCs [41].Furthermore, a positive correlation between seeding density and ALP gene expression has been reported for human umbilical cord MSCs [44], in agreement with the effect observed in this work.Interestingly, hMSCs cultured in osteogenic differentiation medium under shear stress exposure in a 2D system became less responsive in the gene expression of alkaline phosphatase, possibly explaining the dose dependency observed in this study [51].
Col1 was seen not to be affected by the magnitude of flow, while it correlated positively with seeding density.Although collagen matrix synthesis has been reported in numerous studies to be positively affected by flow rate [6,50,52], a possible explanation for the lack of correlation between Col1 and fluid flow magnitude could be the fact that by Day 7, which is the initial time point studied in our work, Col1 expression was already downregulated.Similar Col1 temporal trends parallel to elevated alkaline phosphatase levels have, however, been observed in dynamic culture systems [53].

Highlighting Gene Expression Signature Windows via Response Surface Methodology
The overall aim of this study was to introduce a "windows of operation" methodology as a tool to allow instantaneous visualization of multiple gene expression data related to osteogenic differentiation, thereby aiding decision-making regarding the end-point properties of 3D TE constructs obtained in dynamic culture conditions.By mapping out operating spaces using a DoE approach, user-defined levels of process performance may be finally obtained by linking the effects of process parameters (factors) on pre-selected response variables, in this case, gene expression levels.By using gene specific response surfaces presented in Figure 7 and overlaying them as shown in Figure 8, a case study on the development of a gene expression signature over culture time (days) as a function of flow rate and cell seeding density.By gene expression signature, we refer to a combination of genes that will be expressed by cultured cells to predefined values, where minimum criteria or intermediate ones (between a set minimum and maximum value) may be met.Dashed lines (Figure 8) highlight the values of the two selected process parameters that are required in order to produce TE constructs containing cells that express the selected genes above a desired threshold value.As the culture progresses and the osteogenic differentiation process "matures" from Day 12 to Day 15, it may be seen that the gene expression signature target-values are obtained with a wider combination of culture process parameters, since this is a time-dependent biologic process.This ultimately means that based on the selection of the operating process parameters, TE cultures may be steered towards desired directions where the cell population may possess desirable phenotypic characteristics.Confidence intervals (95%), as seen in Figure 5 could be also incorporated in the outlining of a more restricted, but more reliable window of operation.Recently, a similar approach was used for highlighting "windows of operation" of perfusion bioreactors during early culture time points (up to three days of cell culture), by operating in minimal cell detachment and increased gene expression levels of COX2 and OPN [54].
We do realize, however, that the use of a fractional factorial, rather than full-factorial design, aliases (i.e., does not decouple) the interactive effects of remaining factors to the main effect of a factor not making use of the full extent of the possibilities of such a strategy.However, the systematic investigation of the interplay between culture process parameters and neo-tissue quality attributes, such as cell phenotype, will give insights regarding TE bioprocess development and optimization.Overall, it should be stated that the use of multifactorial DoE approaches and response surface models for process optimization in TE applications has recently gained momentum [12,14,16,40,55].These studies highlight the rising need for improved process understanding, in particular when dealing with inherently highly complex TE systems, meeting, at the same time, regulatory requirements.

Conclusions
In this work, we showed that the use of multifactorial DoE can be employed as a tool to link perfusion bioreactor operating conditions (process parameters, such as seeding density and flow rate) to osteogenic differentiation.Furthermore, by overlaying gene-specific response surfaces, an integrated operating process space was highlighted within which predetermined values of six genes of interest (i.e., gene signature) could be minimally met over the course of a bioreactor culture period.Such an approach will enhance our understanding of the highly complex TE bioprocesses and allow the development of predictive models for in vitro cell behavior and, hopefully, in vivo outcomes.This strategy would thus facilitate the transition from bench to clinic and, ultimately, to the market.

Figure 1 .
Figure 1.The three-factor, two-level fractional factorial design showing the process space investigated and the experimental conditions chosen for this work.

Figure 4 .
Figure 4. Gene expression data for hPDCs cultured in CaP-supplemented medium for all experimental conditions.Data shown are the relative expression of the gene of interest normalized to the HPRT1 housekeeping gene.This data set corresponds to the experimental design employed in this study, as shown in Table 1.Mean ± S.D. error of the mean (n = 3).Unpaired t-test (two-tailed): * p < 0.05.ALP, alkaline phosphatase; Col1, collagen type I; OCN, osteocalcin; OPN, osteopontin; OSX, osterix.

Figure 5 .
Figure 5. Representative images for the midpoint experimental time point (B3) showing collagen matrix content by positive Picrosirius Red staining (A,B) and mineral content by positive Alizarin Red staining (C,D), but also as evidenced by bright field images (E,F).Scale bar: 1 mm.

Figure 6 .
Figure 6.Goodness-of-fit between predicted values obtained via partial least-squares fitting and experimentally measured gene expression values, for the selected genes studied in this study.R 2 values are also given for all plots.

Figure 7 .
Figure 7. Response surface representations showing the effect of seeding density (Seed.Dens.) and flow rate on gene expression levels for all genes studied.Relative expression levels were used to plot the response surfaces shown.The data shown are indicative of Day 14.

Figure 8 .
Figure 8. Development of a selected gene expression signature over culture time (as indicated by the dashed lines) in a perfusion bioreactor setup, as a function of flow rate and cell seeding efficiency.

Table 1 .
(A) DNA content after perfusion bioreactor culture for all experimental conditions; (B) Pareto chart showing levels of significance of process parameters (inputs) on final DNA content (output) of in vitro cultured TE constructs.