Computational model-informed design and bioprinting of cell-patterned constructs for bone tissue engineering

Three-dimensional (3D) bioprinting is a rapidly advancing tissue engineering technology that holds great promise for the regeneration of several tissues, including bone. However, to generate a successful 3D bone tissue engineering construct, additional complexities should be taken into account such as nutrient and oxygen delivery, which is often insufficient after implantation in large bone defects. We propose that a well-designed tissue engineering construct, that is, an implant with a specific spatial pattern of cells in a matrix, will improve the healing outcome. By using a computational model of bone regeneration we show that particular cell patterns in tissue engineering constructs are able to enhance bone regeneration compared to uniform ones. We successfully bioprinted one of the most promising cell-gradient patterns by using cell-laden hydrogels with varying cell densities and observed a high cell viability for three days following the bioprinting process. In summary, we present a novel strategy for the biofabrication of bone tissue engineering constructs by designing cell-gradient patterns based on a computational model of bone regeneration, and successfully bioprinting the chosen design. This integrated approach may increase the success rate of implanted tissue engineering constructs for critical size bone defects and also can find a wider application in the biofabrication of other types of tissue engineering constructs.


Introduction
Bone is one of the most transplanted tissues, with yearly around 2.2 million operations being performed worldwide [1]. Nevertheless, autologous or allogeneic bone grafting procedures have many drawbacks including donor site pain, infections, and disease transmission [2]. Distraction osteogenesis, a bone lengthening procedure, provides a surgical alternative but requires long healing times [3,4]. Therefore, tissue engineering, combining cells with growth factor technology and a biomaterial carrier (scaffold) into living constructs that promote tissue regeneration, is a promising alternative [5]. However, difficulties in healing bone defects persist mainly due to the heterogeneous and often insufficient nutrient and oxygen supply at the site of transplantation leading to cell death and incomplete regeneration/healing [6]. Given the spatial heterogeneity of the healing process, we hypothesize that the tissue engineering construct design must consider spatially heterogeneous (i.e., patterned) initial distributions of biological factors. We therefore propose that a well-designed tissue engineering construct, with an optimized spatial pattern of cell densities in a material matrix, will improve fracture healing in large segmental defects.
However, designing a pattern of living cells which is favorable for a specific biological outcome (e.g., enhanced bone formation) is far from trivial for many reasons including the challenge in understanding the relation between the input pattern of the cells and the biological outcome as well as the many design degrees of freedom possible with, for example, biofabrication technology. Computational modeling provides a tool to qualitatively and quantitatively investigate the relation between the patterned construct and the biological outcome and, based on this knowledge, to design beneficial patterned constructs. Moreover, computational models offer a practical, economical and ethical alternative to in vivo experiments by testing a large number of different conditions and finding out the most promising ones. Many approaches in the literature exist for the optimization of scaffold design parameters related to mechanical and mass transport properties including state-of-the-art finite element analysis and computational fluid dynamical techniques [7]. Byrne et al for example, used a mechanoregulation algorithm based on distortional strain and fluid velocity to model the bone-scaffold system and optimized three design parameters, i.e., scaffold porosity, Young's modulus and dissolution rate, to suit the defect-specific loading requirements [8]. A similar approach was taken by Adachi et al who optimized the scaffold microstructure in order to match the stiffness of healthy bone [9]. However, these computational frameworks for scaffold optimization currently lack a detailed description of the biological aspects of bone regeneration, including spatial patterns of cells and biological signals. Previously, Carlier et al have established a multiscale computational framework of bone regeneration that quantitatively describes the interplay between cells, growth factors, oxygen tension and blood vessels [10][11][12]. In short, this multiscale model combines ten partial differential equations of the taxis-reaction-diffusion type at the tissue level with a discrete agent-based approach at the vascular level, including eight intracellular variables for the endothelial cells. At the tissue level, the model accounts for the various key processes of intramembraneous and endochondral ossification that occur during the soft and hard callus phase of bone fracture healing. At the cellular level, the development of the discrete vascular tree (composed of endothelial cells) is determined by three different processes, i.e., sprouting (the formation of a new branch, headed by a 'tip cell' or 'leader cell'), vascular growth (the extension of the branch due to the growth of 'stalk cells' or 'follower cells') and anastomosis (the fusion of two branches). An anastomosis between blood vessels allows for blood flow and the delivery of oxygen and nutrients. The intracellular level considers a number of molecular players that govern endothelial cell movement. By including a specific description of spatial (bioprinted) cell patterns in the computational model, their effect on the healing of a critical size bone defect can be predicted in silico.
To create cell-patterned structures in the context of tissue engineering and regenerative medicine, bioprinting has taken a revolutionary role in recent years. Bioprinting is a process that allows the spatial positioning of living cells with/without biomaterials (bio-ink) by depositing them layer-by-layer with precision to form tissue and organ analogs [13,14], and its potential has already been demonstrated in several studies such as the biofabrication of aorta [15,16], skin [17], tracheal splint [18], heart tissue [19,20] and bone constructs [21,22]. There are four main approaches [23] to deposit bio-inks: inkjet [24][25][26][27], extrusion [15,28], acoustic [29,30] and laser-based [31][32][33][34]. In inkjet bioprinters, which are analogues to the inkjet printers, the bio-ink is dispensed as drops by the use of thermal or piezoelectric transducers [35]. Microextrusion based bioprinters, on the other hand, extrude bio-ink in the form of continuous stripes rather than liquid droplets [15,16]. Acoustic bioprinters generate acoustic waves to induce the placement of the bio-ink [15,30], and finally laser bioprinters make use of laser beams to transfer the bio-ink onto a substrate [33]. Even though each system has its own advantages and disadvantages, they share certain properties [13,16,36,37]: (1) preprocessing to generate a computer-aided design by employing suitable imaging techniques, (2) printing to produce a physical three-dimensional (3D) construct and (3) post-processing to improve the maturation of the 3D construct where possible. In this work, a commercially available microextrusion bioprinter (Organovo) is modified to dispense cell-laden agarose hydrogels in computationally designed spatial cell patterns. Direct-write bioprinting of cell-laden hydrogels has been previously reported where a mixture of methacrylated ethanolamide gelatin and methacrylated hyaluronic acid or extracellular matrix derived methacrylated gelatin was used [38][39][40]. However, in these studies a photopolymerization step was required prior to printing to obtain a stable construct after printing. In our study, which takes advantage of the thermoresponsive property of agarose, it is possible to obtain structurally stable cell-laden constructs with a high cellular viability. With the modification we introduce to our bioprinting system, it is easy to work with multiple starting materials and thus create gradient constructs.
In this work, we describe a proof-of-concept study that combines computational-informed design with customized bioprinting of cells and biomaterial carriers in the context of bone tissue engineering (figure 1). First, we use the previously established computational model of bone regeneration [10,41] to investigate the influence of continuous stem cell patterns in a bone tissue engineering construct on the bone healing outcome in silico. In a next step, the given constraints on bioprinting resolution are taken into account, and the influence thereof on the bone healing outcome is explored with the computational model. One of the most promising discrete patterns-as identified in silico-is then bioprinted using our bioprinting system. The bio-ink is aspirated from microcentrifuge tubes which can easily be maintained in standard heaters. Hence, the cell-laden hydrogel precursor can be kept at a desired constant temperature throughout the whole preparation process prior to bioprinting which prevents premature gelation as well as exposure of the cells to high temperatures or harmful ultraviolet light. Finally, the bioprinted cellpatterned constructs are characterized to check the designed cell pattern and cell viability.

Computational model
The bone healing outcome of a bone tissue engineering construct with a particular spatial pattern of mesenchymal stem cells (MSCs) was simulated by the previously established multiscale computational model of bone regeneration, which has been validated for both small and large bone defects [10,12]. An overview of the model variables and the initial conditions can be found in figure 2 and table 1.
The geometrical domain of the computational model was deduced from the real callus geometry in a standardized murine fracture model [43]. Due to reasons of symmetry, only one-fourth of the domain was simulated (figures 2(B)-(C)). The ROWMAP solver used to solve the system of partial differential equations requires that the geometrical domain is subdivided into rectangular subdomains, meaning that the domain boundary will be composed of straight lines. As the computational time scales with the number of rectangles, we have previously established a simplified domain geometry that coarsens the smooth outer surface of the periosteal callus into discrete steps, without a loss of accuracy in predicted bone healing outcomes [44]. To simulate the bone regeneration process in a large segmental bone defect, a gap size of 4 mm was used, which is in the same range as other mouse femoral critical defect sizes reported in the literature: 2 mm [45], 3 mm [46], 3.5 mm [47], 4 and 5 mm [43]. The parameter values of the computational model were derived from extensive literature reviews, where possible, or estimated if no relevant data was available [10,41,48,49]. To ensure the existence, uniqueness and non-negativity of the solution, the system of equations was also complemented by suitable initial and boundary conditions (figure 2(C)). MSCs and fibroblasts are released into the callus tissue until post fracture day (PFD) 3 from three possible sources (the periosteum, the surrounding soft tissues and the marrow cavity) [50]. The tip cells are initialized at specific starting positions, reflecting the intact vasculature in the cortical bone and marrow cavity [50]. The osteochondrogenic growth factor was assumed to originate from the fractured bone ends and the cortex away from the fracture site during 5 d and 10 d respectively [51,52]. The complete description of the set of equations, the boundary and initial conditions, the parameter values, implementation details, the underlying assumptions and simplifications can be found in previous publications [10,41,48,49] as well as in the supplementary material.

Implementation of the cell patterns in the computational model
To investigate the influence of a specific bioprinted cell pattern on the bone healing outcome in a large defect size, the initial conditions, describing the pattern of  Table 1. Overview of the model variables and initial conditions. X represents the x-coordinate in figure 2. Peripheral zone ('zone 1B + 3b') and central zone ('zone 1a + 2 + 3a') are explained in the same figure. Bone matrix density 0 mg ml −1 0 mg ml −1 g bc Osteochondrogenic growth factor concentration 100 ng ml −1 10 ng ml −1 g v Vascular growth factor concentration 0 ng ml −1 0 ng ml −1 n Oxygen tension 3.7% + 2.33% (X-2 mm) 3.7% stem cells, were altered. The patterned construct, represented by zones 1a + 2 + 3a of figure 2(C), was modeled as a specific initial spatial pattern of stem cells (figure 3, and table 1). The investigated patterns were continuous or discrete (to account for the bioprinting resolution). The continuous patterns were modeled as a combination of a baseline value and a constant gradient (i.e., cell density is a linear function of a coordinate) in the x, y, or x-y direction (see figure 3). Each grid cell of the domain discretization (see [12] for implementation details) was assigned its local cell density based on its position and the linear function.
Because of the small grid size (25 μm), this resulted in a continuous (smooth) pattern. In order to discriminate between the influence of the pattern and the total average cell density, nine different continuous patterns were calculated for ten different total average cell densities ranging from 10 4 to 10 5 cells ml −1 (figure 3).
As typical bioprinting resolutions are much lower than the applied grid size, such smooth (continuous) patterns cannot be printed. The influence of the bioprinting resolution on the bone healing outcome was therefore explored in silico by coarsening the most optimal pattern, i.e., the continuous pattern of cells was replaced by a discrete pattern consisting of four stripes of 500 μm or eight stripes of 250 μm, thus taking into account the given bioprinting constraints (diameter of the glass capillary) ( figure 3). Here, a unique combination of letters (A-AE) is used to designate specific cases of cell density patterns that are discussed in the following text and figures. A complete list of the linear functions describing the x, y, x-y continuous patterns as well as the values of the individual stripes of the discrete patterns can be found in the supplementary material. In the peripheral callus region adjacent to the construct, i.e., zone 1b and 3b of figure 2(C), the initial conditions were equal to those that led to healing in a small defect (see table 1) [41].

Sensitivity analysis of a discrete cell pattern
In order to explore the large design space of potentially beneficial discrete stem cell patterns, a large sensitivity In order to account for the given printing constraints, the continuous pattern was coarsened to a discrete pattern, consisting of four stripes with a diameter of 500 μm. The bioprinted cell-gradient set is composed of four agarose stripes each containing a different cell density in an ascending order (0.1 × 10 5 , 0.5 × 10 5 , 10 5 and 5 × 10 5 cells ml −1 ). In each petri dish, duplicates of the cell gradient set were bioprinted. A large in silico sensitivity analysis using 'design of experiments' was performed on discrete patterns consisting of eight stripes of 250 μm.
analysis using 'design of experiments' (DOE) was conducted. DOE (or experimental design) is a statistical tool that generates an array of combinations of different parameter values within a predefined parameter space [53]. The sensitivity analysis was done for one discrete pattern that consists of eight 250 μm wide stripes (corresponding to a total width that equals half of the fracture gap) and that are perpendicular to the long bone axis (i.e., pattern in the x-direction; figure 3).
Here, the parameters correspond to the cell densities of the stripes, leading to eight parameters in total (figure 3). The stripe width was chosen to be 250 μm in order to match the bioprinting constraints (equal to the diameter of the small glass capillary of the bioprinter; see section 2.6) and in order to have more design freedom (compared to the 500 μm glass capillary, which would result in only four stripes and therefore parameters). Next, the computer simulations were run for 100 combinations of these parameter values, and the results were analyzed. In uniform designs, the parameter combinations are spread out over the parameter space as uniformly as possible. They are found to be efficient and robust but computationally very demanding [54]. Since uniform designs cope well with the addition and removal of parameter combinations, a uniform design with 100 runs was chosen for the sensitivity analysis. The parameter space was determined by varying the stem cell densities of the eight stripes between 0 and 1 × 10 5 cells ml −1 . The uniform designs were generated and analyzed with JMP (SAS Institute Inc.). The cell center of mass (i.e., center of mass of the cell density distribution) of each pattern was calculated as follows: with M: the sum of the densities of all stripes, r : i the cell density of an individual stripe, x : i the x-coordinate of the midline of an individual stripe and n: the number of stripes.

Evaluation of the computational model results
In this study, we focused on the predicted bone tissue fraction after 90 d as the main model outcome. Tissue fractions were calculated by first binarizing the spatial predictions using tissue-specific thresholds (where 0 means that no tissue is present while 1 means that the tissue is present in a grid cell). Subsequently, an equal weight is assigned to the different tissues, i.e., fibrous matrix density (m f ), cartilaginous matrix density (m c ) and bone matrix density (m b ). If a grid cell contains multiple tissues, the area of that grid cell is divided by the number of tissues in the final calculations of the tissue (area) fractions. The time to healing is defined as the earliest time point at which the bone tissue fraction reaches its maximal value (and steady state). When the bone tissue crosses the fracture gap, a bone bridge, connecting the original bone ends, is established, and this is defined as reaching a clinical union. The time efficiency of a particular pattern is defined as the ratio of the maximal bone tissue fraction to the time to healing. The model outcomes (i.e., bone tissue fraction and clinical union) of all continuous and discrete patterns can be found in the supplementary material.
2.5. Cell culture NIH 3T3 murine embryonic fibroblasts (ATCC) were selected for this proof-of-concept study because of their ease of use. They were grown in Dulbecco's modified eagle medium (DMEM; Sigma) supplemented with 10% fetal bovine serum (FBS; Sigma) and 1% penicillin-streptomycin (Gibco). Cells were passaged twice a week and maintained in a humidified incubator at 5% CO 2 and 37°C.

3D bioprinting system and bio-ink preparation
In this work a modified bioprinter (Organovo) was used. It was operated inside of a standard biosafety cabinet for sterile working conditions [16]. This bioprinter has two deposition heads, one of which can print hydrogels while the other one can print cells ( figure 1(B)). The deposition heads can be equipped with 250 or 500 μm in diameter and 85 mm long glass capillaries which dispense hydrogels/cells as cylindrical stripes. In the original design of the system, the glass capillaries aspirate hydrogel from a glass vial placed in the gel chamber ( figure 1(B)) that can be heated to a desired temperature. The usage of glass vials causes several problems especially in the case of handling thermoresponsive hydrogels. The glass vials do not fit in commercially available heating/cooling devices, and if multiple cell/material formulations are to be printed within the same construct, it is difficult to preserve the desired viscosity of the hydrogels. Here, we modified the system and adapted it to use microcentrifuge tubes by introducing physical changes to the bioprinter platform and changing the script provided in the original bioprinter software. This modification allowed us to aspirate materials from microcentrifuge tubes instead of glass vials. We used a heater shaker to prepare and also to maintain the bioink to be printed at the desired temperature and rotation speed ( figure 1(B)). In this way, it was possible to prepare as well as to maintain multiple samples, with varying properties, at a constant temperature thereby at the desired viscosity.
To prepare the bio-ink, in other words, the cellladen hydrogel, the following steps were performed. First, agarose hydrogels were prepared by mixing 3 g of low gelling temperature agarose (Biozym Sieve 3:1), which forms an irreversible hydrogel after cooling below 38°C, with 100 ml 1X phosphate buffered saline (PBS; Hyclone by Thermo Scientific). The resulting mixture was autoclaved under standard liquid sterilization conditions and kept in liquid state until usage. Next, cell cultures were washed with 1X PBS, trypsinized for 5 min and centrifuged at 1100 rpm for 5 min. The resulting cell pellet was resuspended in growth medium, stained with trypan blue (Sigma) and counted with a hemocytometer. Third, previously prepared 3% liquid agarose was transferred into microcentrifuge tubes and kept on the heater shaker at 38°C, 300 rpm until its initial temperature (∼60°C) dropped to 38°C at which it was still preserved in the liquid state. Finally, the required numbers of cells were mixed with the liquid agarose gel precursor. Therefore, a uniform cell distribution was obtained before the agarose solidified at ∼37°C, and in addition, cells were not exposed to high temperatures. The following four different cell concentrations were used to prepare the cell-laden agarose bio-ink based on the outcome of the computational model: 0.1 × 10 5 , 0.5 × 10 5 , 10 5 and 5 × 10 5 cells ml −1 .

Bioprinting of cell-laden agarose hydrogels
We chose the 500 μm glass capillaries to print and characterize a specific cell-gradient construct, predicted by the computational model to result in a union. Hence, 500 μm diameter sized dispensing glass capillaries were used to aspirate the cell-laden agarose from the gel chamber which was adapted to hold microcentrifuge tubes. The capillary loaded with cellladen agarose mixture was cooled to 4°C by being submerged in PBS, and the cell-laden agarose hydrogel was then mechanically extruded directly on glass bottom petri dishes. The resulting cell-laden agarose stripes preserved their cylindrical shape and were each 1 cm long and 500 μm in diameter ( figure 3(B)). They contained different cell concentrations and were located next to each other so that a gradient was obtained ( figure 3(B)). Each petri dish contained duplicates of the gradient set which was composed of four cell-laden agarose stripes with distinct cell concentrations. The bioprinted constructs were then incubated for up to 3 days submerged in growth medium at 37°C in a humidified atmosphere containing 5% CO 2 . The bioprinted construct is thermally irreversible and the shape of the construct does not change in the incubator. The cell-laden agarose stripes were analyzed immediately after bioprinting (Day 0) as well as 1 and 3 days of culture (Day 1 and Day 3, respectively).

Characterization of the bioprinted cell-gradient constructs
Prior to bioprinting, the viability of the cell suspension at Day 0 was evaluated by performing a trypan blue exclusion test. The cell suspension was mixed with trypan blue, and the cells were visually examined and scored whether they took up (dead) or excluded the dye (live). The fraction of viability values was obtained by dividing the number of trypan blue negative cells by the total number of cells. Cell viability within the bioprinted cell-gradient constructs was assessed by applying a live/dead fluorescence assay at Day 0, Day 1 and Day 3. The constructs were stained with 2 μM calcein-AM (Invitrogen, green fluorescence) for 20 min, briefly washed with 1X PBS and then incubated in 0.5 μM propidium iodide (Invitrogen, red fluorescence) for 5 min. They were stained and imaged directly in glass bottom petri dishes without being disturbed. For imaging an inverted confocal microscope (Zeiss LSM710) was used. Tiled z stacks were performed which resulted in images of cylindrical cell-laden agarose stripes with a diameter of 500 μm and a length of 5.65 mm. This cylindrical volume was used as the unit volume for the quantitative analysis. The 3D reconstructed images were quantitatively analyzed by Image J 1.48v software [55]. The 3D object counter plug-in of Image J software was employed to score the number of green/red fluorescent cells encapsulated in the unit volume of agarose hydrogel [56]. The fraction of viability per each cellladen agarose stripe with varying cell density was calculated by counting the number of live cells (green) as well as dead cells (red) and dividing the number of live cells by the number of total cells (live + dead). The cell viability of the whole construct, consisting of 4 cell-laden agarose stripes with 0.1 × 10 5 , 0.5 × 10 5 , 10 5 and 5 × 10 5 cells ml −1 , was calculated by summing up the viability values of these 4 cell-laden stripes and dividing the resulting number by 4.

Statistical analysis
All data is presented as the mean ± standard error of the mean. The mean values represent the average values of three independent experiments which were done in duplicate. For statistical analysis one/two-way analysis of variance followed by a Tukey post hoc test and a student t-test were performed. Statistical significance is defined as * p < 0.05, and ** p < 0.01.

Results
3.1. Patterned constructs are able to enhance bone regeneration in silico Figure 4 plots the bone tissue fraction at post fracture day 90 (PFD 90) as a function of the total average cell density for eight different linear stem cell distributions and a uniform pattern. The particular pattern-total average cell density combinations that resulted in a union are indicated with an arrow head (see figures 4 and 6). Figure 4 shows a highly nonlinear relationship between the bone tissue fraction at day 90 and the total average cell density. Both at low (0.1 × 10 5 cells ml −1 ) and high (1 × 10 5 cells ml −1 ) total average cell densities the bone healing was impaired, indicated by the low bone tissue fraction, regardless of the spatial stem cell pattern. Only total average cell densities between 0.3 × 10 5 and 0.8 × 10 5 cells ml −1 resulted in complete healing for some of the stem cell patterns (indicated by the arrow heads). None of the cell densities led to union in the case of a uniform distribution.
At the intermediate cell densities, the bone tissue fraction was substantially influenced by the initial stem cell pattern. Pattern 2 (indicated by the blue line in figure 4: linear distribution with highest cell density near the bone ends) is the most robust with respect to changes of the total average cell density since it provides successful healing for four different total average cell densities (i.e., 0.5, 0.6, 0.7 and 0.8 × 10 5 cells ml −1 ). Interestingly, successful healing (indicated by the arrow heads in figure 4) can be obtained for a variety of patterns, but for a given pattern the amount of bone formation strongly depends on the total average seeded cell density (see, e.g., pattern 3). For every pattern, except for pattern 8, there is at least one pattern-cell density combination that led to a complete healing of the bone defect. The uniform distribution of cells (indicated by the red, dotted line of figure 4) neither reached high bone tissue fractions nor resulted in a clinical union, independent of the total average cell density. However, at particular total average cell densities (e.g., 0.6 × 10 5 cells ml −1 ) the uniform distribution did perform better than some stem cell patterns (e.g., patterns 3, 4, 5, 8 and 9). It is also interesting to note that a high bone tissue fraction does not necessarily imply successful healing since the bone tissue fraction does not account for the spatial distribution of the bone and, as such, the bridging of the bone defect. The spatial distribution of the bone formation at PFD 90 is shown for selected cases in figures 5 and 6. An example of this phenomenon is given in figure 5 for pattern 7, which resulted in a union at 0.7 × 10 5 cells ml −1 (case L, 85% bone) but a non-union at 0.4 × 10 5 cells ml −1 (case N, 82% bone) even though the bone tissue fractions were very similar. Similarly, pattern 5 resulted in a union at 0.5 × 10 5 cells ml −1 (80%) whereas pattern 7 resulted in a non-union at 0.4 × 10 5 cells ml −1 (82%) although the latter had a higher bone tissue fraction at PFD 90 ( figure 4).
A more detailed look at the specific pattern-total average cell density combinations that resulted in a clinical union is given in figure 6. The time required for  . The predicted amount of bone formation at day 90 (× 0.1 g ml −1 ) for a specific stem cell distribution at day 0 (× 10 5 cells ml −1 ) in a large segmental defect. Note that the amount of bone formed during healing was similar in case L and N but that the spatial distribution of the formed bone was different, leading to bone bridging (and consequently a clinical union) in case L only.
healing varied between 75 and 90 days where a higher bone tissue fraction often corresponded to a longer healing time. In this light, cases A and J were the most efficient as they reached a bone tissue fraction of 100% after 83 days (figure 6). Note again the importance of the bone formation spatial distribution in assessing the potential of a specific stem cell pattern. For example, even though case D had a higher bone tissue fraction than case L, according to the computational model, very little bone tissue crossed the fracture gap (limited bone bridging) in case D, possibly resulting in a reduced mechanical strength of the bone. Consequently, the clinical outcome of case L could be much better than that of case D, especially when the reduced time to healing and the increased efficiency are also taken into account. Figure 6. The predicted spatial distribution of bone formation at day 90 (× 0.1 g ml −1 ) for the specific pattern-total average cell density combinations at day 0 (× 10 5 cells ml −1 ), that resulted in a clinical union of the large segmental defect in silico (see also figure 4, arrow heads). The pattern number of the nine patterns is explained in figure 4. A complete list of the equations describing the x, y, x-y continuous patterns can be found in the supplementary material.

Influence of printing resolution on successful computational cell-gradient patterns
The results above revealed that particular patterned constructs enhanced bone regeneration. These simulated patterns (i.e., initial cell densities) were very smooth (continuous) and did not take into account any constraints given by the bioprinting resolution (250 or 500 μm in our case). Interestingly, when a pattern was coarsened, in order to account for this constraint, the bone healing outcome was highly affected (figure 7). As an example, lowering the resolution of the stem cell pattern from 25 μm (i.e., the spatial resolution given by the grid size) to 250 μm (eight discrete stripes), resulted in less bone formation for case O with respect to case B. On the other hand, a further reduction in the resolution (four stripes of 500 μm wide) increased the amount of bone formation (figure 7, case P) and even resulted in a complete healing at a specific combination of four cell densities (figure 7, case Q). In contrast, the lower resolution of case R yielded more bone formation than the resolution of case F. A further reduction in the resolution decreased the amount of bone formation (figure 7, case S) or resulted in a non-union ( figure 7, case T). According to the model, case Q, consisting of four stripes containing 0.1 × 10 5 , 0.5 × 10 5 , 1 × 10 5 and 5 × 10 5 cells ml −1 , resulted in complete healing of the large segmental defect; therefore, a micro-extrusion based bioprinter was used to bioprint this promising discrete pattern and characterize it in vitro (see section 3.4).

Computational models can aid the design of novel patterned constructs
Inspired by the results of figure 7, a large sensitivity analysis using DOE was performed to further investigate the influence of eight individual stripe cell densities on the bone healing outcome. Figure 8 shows the bone tissue fraction at PFD 90 as a function of the cell center of mass (A) and the total average cell density (B). First, it can be noted that the simulation results span a wide range of bone tissue fractions (50%-100%). Second, only a small number (11 out of 100) of stem cell patterns lead to a clinical union (indicated by the dots within the dashed box). Moreover, the successful patterns (indicated by the dots within the dashed box) tended to have a total average cell density in the lower half of the explored density range and a center of mass closer to the edge of the defect. An overview of the discrete patterns that resulted in a clinical union is given in figure 9. The time to healing varied between 77 and 90 days. In this light, discrete case AE was the most efficient as it reached a bone tissue fraction of 100% after 77 days. It is also interesting to note that this case was more efficient than any of the continuous patterns previously investigated ( figure 6).

Characterization of 3D patterned cell-gradient constructs
Confocal fluorescent images demonstrated a high cellular viability among all samples at all time points (figures 10(A) and (B)). In addition, the preservation of the cylindrical architecture of the cell-laden agarose as well as the presence of a cell-gradient were observed. The number of live (green) and dead (red) cells were then counted in each unit volume of cell-laden agarose stripe seen in figure 10(A). High cell viability was observed at all time points. For the whole construct there were no statistically significant differences between the viabilities at all time points. At Day 1 and Day 3, the majority of the cells were still viable resulting in 93% viability while at Day 0 the viability was 90% ( figure 10(C)). Regarding the individual viabilities of each stripe with varying cell concentrations, values ranging between 85% and 97% were obtained (figure 10(C)). A statistical analysis indicates that cell density has a significant effect on cell viability (figure 10(C)). At early timepoints cell-laden agarose stripes with lower cell densities have lower cell viabilities.
We also quantitatively analyzed the presence of the gradient composed of four different cell concentrations. Figure 11 compares the measured total cell numbers with the expected ones. The expected values were calculated based on the volume of the cylindrical agarose stripe and the concentration of the cell suspension used to prepare the cell-laden agarose. Our results show that the observed values were very close to the expected ones, and no significant differences were found between observed and expected data for any of the cell densities ( figure 11).

Discussion
In this study, we provide a novel strategy for bone tissue engineering in the context of large bone defects and demonstrate the feasibility of bioprinting constructs with computationally designed spatial patterns of cells. We used a well-established computational model of bone regeneration [10,41] to show that, compared to uniform constructs, spatially patterned constructs are able to enhance bone regeneration in silico. Interestingly, we also found that discretization of these patterns, which is required to design bioprintable constructs, has an important effect on the bone regeneration outcome. Further, we successfully obtained cell-gradient constructs with high cellular viability by using a modified bioprinting system.
An important direction within the (bone) tissue engineering field has been to strive towards a uniform cell distribution within a tissue engineering construct before implantation [57,58]. It is considered that a uniform cell distribution would promote uniform extracellular matrix distribution and thus result in uniform tissue formation and integration with the native tissue [6]. However, recent advances in the field have been focusing on the creation of discrete or continuous gradients of biological molecules and/or cells in order to mimic their spatiotemporal variation during various healing processes (e.g., neovascularization [59]) or to recapitulate the natural gradients present under physiological conditions (e.g., the osteochondral interface [60]). Examples of tissue engineering constructs presenting a gradient in pore size [61], mechanical properties [62,63] or biochemical composition [64] as well as gradients in growth factor concentrations [65,66] have been reported in literature. Similarly, in this study, we successfully designed, bioprinted and characterized a cell-gradient construct. Our in silico predictions indicate that the spatial heterogeneity of the natural healing response with, e.g., more oxygen near the existing and developing blood vessels than in the center of the defect (see supplementary figure S1), requires a spatially heterogeneous (i.e., patterned) initial distribution of stem cells. Indeed, the computational results clearly demonstrate that patterned constructs are able to enhance the bone regeneration compared to uniform constructs. Furthermore, the regenerative potential of the bone construct was found to be dependent on the total average cell density, the specific stem cell pattern as well as the bioprinting constraints such as the bioprinting resolution. In particular, only total average cell densities between 0.3 × 10 5 and 0.8 × 10 5 cells ml −1 resulted in a union for some of the continuous stem cell patterns (indicated by the arrow heads, figure 4). This result corresponds to a previous finding [11] and can be explained as follows. On the one hand, a low initial total average cell density reduces the biological potential of the fracture site as fewer cells can contribute to the bone regeneration process. On the other hand, a high initial total average cell density will worsen the detrimental hypoxic conditions in the center of the defect due to an increased amount of oxygen consumption in the absence of new blood vessels that start developing at the edge of the bone defect (see also discussion below and supplementary figure S1). At the intermediate cell densities, the bone tissue fraction was substantially influenced by the initial stem cell pattern, with pattern 2 being the most robust with respect to changes of the total average cell density. Similar to the reasoning regarding the low and high total average cell densities above, the vascularization of the host environment will influence the success of particular gradient constructs. More specifically, pattern 2 is characterized by a high cell density near the blood vessels of the surrounding soft tissues, ensuring adequate oxygen delivery and survival of the implanted cells. Moreover, in the central defect region, the hypoxic conditions are limited due to the low cell density, ensuring cartilage formation and adequate endochondral ossification. Interestingly, a particular discretization of pattern 2 (case O), resulted in a nonunion due to hypoxia and cell death (see supplementary figure S2) whereas other, less intuitive stripe combinations, not resembling the continuous spatial distribution of pattern 2 did result in complete bridging (see figure 9  and discussion below). Moreover, the main direction of bone formation is axial (from right to left, see supplementary figure S3) independent from the initial cell gradient. The particular distributions of bone density (see figures 5 and 6) arise due to the heterogeneous distributions of cartilage, which are in turn influenced by the local survival of the chondrocytes arising from the initial stem cell pattern and oxygen distribution.
In order to study the influence of the bioprinting resolution on the bone healing outcome in more detail, the continuous (smooth) cell pattern was replaced by discrete patterns of either 500 μm (four stripes) or 250 μm (eight stripes) resolution. The results demonstrate that successful patterns can be designed at both printing resolutions albeit more flexibility is given by the eight stripe-pattern (figures 7 and 9). Therefore, the large design space of the eight-stripe cell pattern was explored in more detail using DOE. Only 11% of the 100 different stripe density combinations resulted in a clinical union. Interestingly, these successful patterns tended to have a total average cell density in the lower half of the explored density range and a center of mass closer to the edge of the defect. In the in silico model the blood vessels, growing from the surrounding soft tissues in the large bone defect, are initialized at the bone ends (see supplementary figure  S1). As such, the new vasculature grows axially, from the edge of the large segmental defect towards the center. If we contrast the successful patterns with the vascularization of the host environment, we might conclude that the cell density of successful patterns is in general higher in regions that are close to blood vessels in the surrounding soft tissues. In other words, the successful patterns might satisfy specific design guidelines such as a cell center of mass as calculated by the current model specifications > 800 μm (resulting in a higher cell density near blood vessels), although further research is warranted to specify these particular design guidelines and their respective thresholds. Taking all the in silico results together, it appears that the computational model predicts the complete healing of a large segmental bone defect for multiple, very distinct cell patterns at different bioprinting resolutions (figures 6, 7 and 9). Moreover, the model is able to identify less intuitive patterns that perform better than the more obvious ones (compare cases Z and AE in figure 9). We may therefore state that the proposed approach can aid in the design of beneficial cell patterns for bone tissue engineering constructs.
After determining a favorable, discrete pattern in silico, we successfully bioprinted constructs with the designed cell-gradient pattern. The concept of pattern printing has been demonstrated by earlier studies [67][68][69]. For example, in one study collagen was printed in patterns by using an inkjet printer [67] while in another study serial dilutions of a protein of interest encapsulated in a hydrogel were 3D printed [68]. In these reports the main objective was to study cell behavior in a gradient environment. Similar approaches were also used to direct stem cell fate by patterning biomolecules like growth factors [36]. In this work, we bioprinted constructs with a gradient in cell density, which to our knowledge has not been reported previously. In addition, the successfully bioprinted gradient construct is composed of low cell densities. Such an approach might especially be useful for studies where precious cell sources like stem cells, available only in limited numbers, are to be used. The use of cell-laden hydrogels has been explored in earlier studies where either direct printing or formation of 3D tissue models by molding was performed [39,40,70,71]. In these studies the cell densities used were 5 × 10 5 , 10 6 and 2 × 10 6 cells ml −1 respectively while we attained 0.1 × 10 5 cells ml −1 in our lowest cell density stripe. Our quantitative analysis, indeed, shows that the obtained cell-gradient values were very close to the expected ones. Although it is not possible to obtain cell densities with the precision of a single cell in cell-laden hydrogels with the current technology, we have demonstrated for the first time the feasibility of bioprinting both very low and high densities of cells encapsulated in a hydrogel within the same construct.
Preserving cell viability is one of the major concerns related to bioprinting. Our results demonstrated a high cellular viability over a culture period of 3 days, with values reaching 90% or more over the entire construct. These values are similar to previous reports about NIH 3T3 cell encapsulation in non-bioprinted agarose hydrogels [47]. It should be noted that in our study high cell viability values could be maintained despite the fact that four different hydrogel solutions with varying cell concentrations were used within the same construct. This could be attributed to the modification we applied to our bioprinting system that allowed us to handle four different cell-laden hydrogels at temperatures close to normal body temperature. With this developed method, it is also possible to combine different cell populations within the hydrogels and bioprint them within the same construct by simply alternating the bio-ink vial (microcentrifuge tube in the gel chamber). Interestingly, our data also showed that there is a significant effect of cell density on cell viability at early timepoints. In particular, in cell-laden agarose stripes with lower cell densities, lower cell viabilities were observed at early timepoints. The sparse distribution of cells might have a role on cell viability which can be explored further in future studies. One of the limitations with the presented approach is that the cell viability drops considerably if the same cell-laden agarose solution is used for bioprinting multiple samples (data not shown). This is mainly because cells remain suspended in agarose without access to culture medium for long periods of time. In order to overcome this problem, we prepared fresh cell-laden agarose solutions immediately before the printing process for every petri dish to decrease the length of time that the cells spent in the solution. We obtained a noticeable improvement in viability in this way.
The designed constructs were printed as a proof of concept that the gradient structures coming out of the model results could be printed using capillary-based extrusion. The designed gradient structures can be bioprinted in a layer-by-layer fashion where each layer can have various (2D) designed gradient structures resulting in a complex 3D gradient structure (see supplementary figure S4). The height of the constructs was limited for practical reasons. Larger constructs would not be able to be imaged using confocal microscopy due to the limited imaging depth, and further, a necrotic core may be formed in vitro due to limited diffusive transport in hydrogels. Further upscaling to centimeter scale defects in large animals or humans is likely to require additional modification to the design strategy, to avoid mechanical collapse and to overcome mass transport limitations of relatively soft, limitations of relatively soft hydrogels with low permeability. Recent work has demonstrated the potential of hybrid printing strategies, combining stiff polymers like polycaprolactone with cell-laden hydrogels and at the same time integrating microchannels for improved mass transport. Mechanical stability and tissue formation has been shown for hybrid construct sizes up to one centimeter [72]. In the future, such hybrid printing strategies could be combined with computational modeling to further optimize construct design.
By combining the in silico and in vitro results, we can state that the proposed integrative approach allows us to identify multiple cell patterns that result in a union as well as to successfully bioprint them. A next important step towards clinical application of patterned, bioprinted tissue engineering constructs is the in vivo validation of the in silico predictions, in order to find out whether the identified in silico patterns truly perform better compared to uniform constructs in an in vivo setting (e.g., in a mouse model). Moreover, in vivo experimentation could potentially aid in discriminating within the large set of successful, albeit very different and non-intuitive, patterns by defining additional outcomes such as mechanical functionality. In the future, the predicted bone distributions could serve as input for a finite element analysis to predict the corresponding stiffness at the time of bridging, offering an additional measure of the quality of bridging. Together with the in vivo validation, this could enhance the selection of beneficial cell patterns. Finally, in case of discrepancies between the predicted and measured in vivo outcomes, the in vivo validation would allow for the identification of any shortcomings of the proposed approach, e.g., related to simplifications that are inherent to any computational model (see also [10][11][12] for a more elaborate discussion on model simplifications).

Conclusion
In this study, we describe a novel integrated approach for the biofabrication of bone tissue engineering constructs by combining a computational model of bone regeneration with bioprinting technology. Our in silico results indicate that spatially patterned constructs can enhance bone regeneration compared to constructs with a uniform cell distribution. In addition, the discretization of the pattern, which is required for the bioprinting process, has a substantial effect on bone regeneration and should therefore be accounted for during the design process. Following the model-informed design, we successfully bioprinted patterned constructs and demonstrated that the gradient pattern and cell viability were maintained for 3 days of culturing. This integrated approach may ultimately contribute to the treatment of critical size bone defects and can also find a wider application in the biofabrication of other types of tissue engineering constructs.