Simulations of 3D bioprinting: predicting bioprintability of nanofibrillar inks

3D bioprinting with cell containing bioinks show great promise in the biofabrication of patient specific tissue constructs. To fulfil the multiple requirements of a bioink, a wide range of materials and bioink composition are being developed and evaluated with regard to cell viability, mechanical performance and printability. It is essential that the printability and printing fidelity is not neglected since failure in printing the targeted architecture may be catastrophic for the survival of the cells and consequently the function of the printed tissue. However, experimental evaluation of bioinks printability is time-consuming and must be kept at a minimum, especially when 3D bioprinting with cells that are valuable and costly. This paper demonstrates how experimental evaluation could be complemented with computer based simulations to evaluate newly developed bioinks. Here, a computational fluid dynamics simulation tool was used to study the influence of different printing parameters and evaluate the predictability of the printing process. Based on data from oscillation frequency measurements of the evaluated bioinks, a full stress rheology model was used, where the viscoelastic behaviour of the material was captured. Simulation of the 3D bioprinting process is a powerful tool and will help in reducing the time and cost in the development and evaluation of bioinks. Moreover, it gives the opportunity to isolate parameters such as printing speed, nozzle height, flow rate and printing path to study their influence on the printing fidelity and the viscoelastic stresses within the bioink. The ability to study these features more extensively by simulating the printing process will result in a better understanding of what influences the viability of cells in 3D bioprinted tissue constructs.


Introduction
3D bioprinting has the potential to revolutionise medicine by enabling biofabrication of patient specific human tissues, where cell-laden bioinks are deposited with controlled spatial distribution. The interest in 3D bioprinting has led to a vast development of biomaterials which can be printed with living cells, so called bioinks [1][2][3]. The required properties of a bioink are highly governed by the targeted tissue since it defines the selection of cells, mechanical properties and architectural design. In addition to acting as an ideal scaffold for cells, bioinks must also be printable [4]. In contrast to seeding cells on 3D scaffolds, bioinks can be combined with cells and printed in one step. Several 3D bioprinting techniques have been developed such as ink-jet printing [5], laser-induced forward transfer [6,7], microvalve-based [8] and extrusion-based bioprinting [9,10]. Extrusion-based bioprinting can also use viscous bioinks, which makes it the most versatile with regard to types of printable bioinks [11,12]. Hydrogels are preferred as a bioink since they provide, in similarity to the extracellular matrix of native tissues, an aqueous environment and structural support, which is beneficial for the cell survival and proliferation [13]. A wide range of hydrogels have been evaluated as bioinks; collagen [14,15], alginate [16], gelatin [17,18], gelatin-methacrylate [19,20], nanocellulose [21,22], fibrin [23,24] and hyaluronic acid [25]. Within each category of hydrogels there is a vast amount of bioink compositions where hydrogels have been modified or combined with other materials to increase either their printability or cell interactions. An example of this is alginate which has been used at concentrations ranging from 0.8 to 10 wt%. It has been modified with RGD peptides [26] and furthermore has been combined with nanocellulose to increase the viscosity [21,27] or with gelatin-methacrylate for UVcrosslinking. In order to fulfil the multiple requirements of bioinks, scientists are developing new materials and bioink compositions. In the majority of studies that involve development of new bioinks, the cell viability, proliferation, and growth of cells in the bioink is evaluated thoroughly while the assessment of printability is left at verifying that the ink can be deposited by a bioprinter and that a gridded scaffold can be printed. This is partly due to a lack of standardised methods for evaluating printability of bioinks which prevents the comparison of different bioink compositions. Bioinks are most commonly evaluated experimentally which includes testing a series of printing parameters (flow rate, printing speed, layer height) for various print designs such as printing of lines, grids or cylinders [24,26,28,29]. The printability is then determined with regard to the ratio of line width to nozzle diameter, the amount of layers until collapse, or the curvature of printed lines. Determining the printability by printing trials is timeconsuming as well as material consuming. For a bioink that contains cells, the prepared ink is valuable and testing of the bioinks printability must be kept at a minimum. Moreover, after mixing cells with bioinks there is no time for extensive evaluation because the cell viability is at risk [28,30,31]. This leads to many bioinks containing cells being printed at printing parameters that have not been optimised for that specific bioink composition and therefore the printing fidelity is sacrificed. Consequently, the desired architecture is not achieved which may be catastrophic because the printed tissue is not able to function properly [32]. The printed lines may be too thin causing the structures to break, or too thick, hindering nutrients and oxygen to reach the cells in the bioink. Printing with bioinks containing cells also complicates the evaluation of printability since the viscosity of the bioink may be lowered if diluted by the addition of cell suspension. In conclusion, there is a challenge in determining whether a difference in printing fidelity depends on the material itself, on the printing parameters, or on the addition of cells. Moreover, the printability could be impaired by clogging due to either inhomogenous inks, uneven flows or the design and material of the nozzle tips [1,33,34]. As a complement to experimental testing, an increasing amount of studies are relying on modelling and simulation as part of the 3D bioprinting process and evaluation of bioinks. Most of these studies focus on predicting the shear forces subjected to the bioink in the printer nozzle [27,[35][36][37]. Studying the forces in the printer nozzle has been essential for understanding whether the cells may be subjected to forces harmful for the cell viability. Attempts on predicting the mechanical stability of the final construct, based on material properties and the computer aided design have also been presented [1,38]. These simulations assume that the printing process results in a printed construct identical to the used computer aided design, which is seldom the case when printing with hydrogels in relation to printing of solid materials such as metals or plastics. By simulating the dispensing step of the printing process, the properties that influence specific features of the printing process could be identified and would limit the time and cost of bioink evaluation. Likewise, simulation of the actual printing process, outside of the printer nozzle, is expected to be important since cells used in the bioink may be subjected to shear forces and viscoelastic forces during the deposition of bioink in addition to forces that cells are subjected to inside the nozzle. As described, 3D bioprinting is a process with several parameters that influence the final result and therefore requires a large amount of resources to optimise the printing process experimentally. To reduce the cost and time consumed in experimental evaluation, a complement would be to use computer based simulations. Thereby the printability could be predicted and the simulation could aid in determining the optimal printing parameters prior to printing experimentally. This paper presents how a CFD simulation tool, IPS IBOFlow, can be used to predict the printing process of hydrogels used for bioinks. IPS IBOFlow, which utilises a viscoelastic rheology model and a novel surface tension model, simulates the deposition of bioink and the final shape of the printed material. This work focuses on verifying the simulation tool and showing how simulation could be utilised for evaluating the printing process of bioinks. For instance, by simulating the printing process, parameters that influence the printability can be isolated. Thus, the viscoelastic behaviour of a bioink could be studied separately from the printing speed or the nozzle height, to evaluate what influences the printing resolution and the viscoelastic stresses.

Materials
The bioinks presented and compared in this paper are based on cellulose nanofibril (CNF) dispersions in water and exhibit distinct differences in flow behaviour. The first bioink, ink 6040, is based on CNF dispersions and alginate, and is a less viscous ink with characteristics governed by surface tension forces. CNF was kindly provided by RISE Bioeconomy (Sweden) and alginate (SLG100, Mw=150-250kDA, above 60% α−1− glucuronic acid ) was purchased from FMC Biopolymers (Norway). 3 wt% CNF was mixed with 3 wt% alginate at a dry weight ratio of 60:40. Ink 6040 has been studied within the tissue engineering field and shown promising results as an ink with regard to printing and cell survival [21,22,39]. During processing of CNF from pulp, the pulp is subjected to either enzymatic or chemical pretreatment to increase the efficiency of the mechanical fibrillation and high pressure homogenisation. Ink 6040 was prepared from enzymatically pretreated CNF [40]. The second bioink, 4% CNF ink, was used as a comparison to ink 6040 due to its difference in viscoelastic properties. It consisted of a 4 wt% CNF dispersion in water, and is a more viscous bioink with characteristics driven by viscoelastic stresses. Instead of being enzymatically pretreated, the 4 wt% CNF had been pretreated by carboxymethylation [41] which leaves the CNF more charged, and thereby less aggregated and more dispersed in water.

Rheology measurements
Amplitude sweeps and oscillation frequency measurements were conducted for all inks using a Discovery Hybrid Rheometer (TA instruments, UK) with a 20mm plate-plate geometry at a gap distance of m 250 m. Amplitude sweeps were measured at 1 Hz and controlled stress from 1 to 100 Pa. The frequency sweeps, shown in figure 2, were measured within the linear viscoelastic region from 0.01 Hz to 100 Hz.

Experimental setup
All inks were printed with the same pattern (figure 1), nozzle geometry, height of nozzle and printing speed.
Tissue engineering scaffolds with a porous architecture are preferred to allow for nutrients and oxygen to reach the cells embedded in the bioink [42]. Therefore, 3D bioprinting of grid-like structures are 3D bioprinted when assessing bioinks and their compatibility with cells [43][44][45]. To capture the most important features of 3D bioprinted grids, the pattern was designed to include straight paths, corners and cross overs. A conical nozzle with an outlet diameter of 0.42 mm was used. The distance between the outlet of the nozzle and the printing plate was 0.4 mm for the first layer (continuous line in figure 1) and 0.8 mm for the second layer (dotted line in figure 1). The printing speed, which is the movement of the printer head in the X-and Y-direction, was set to -10 mm s 1 . The flow rate is controlled by adjusting the applied pressure. The pressure was chosen so that a continuous flow was obtained and so that all inks were printed with a flow rate of 90-100 mg min 1 . This applied pressure varied for each ink due to the difference in viscosity and yield stress of the inks. Generally, an ink with a higher yield stress requires a higher pressure force to overcome the fluid stresses and accelerate the fluid ink. The specific flow rate for each ink was determined by weighing the mass of ink dispensed during 30 s at the pressure required to reach the aimed flow rate. The applied pressure and corresponding flow rate (averaged from four measurements) for each ink are presented in table 1. Directly after printing, the printed grid was photographed for comparison with the results from the corresponding simulation. Images of the printed grid from the top and from the side were used to measure the width (n=50) and height (n=20) of the printed lines by ImageJ software [46].

Numerical methods
The simulation software used is IPS IBOFlow, which is an incompressible, segregated Navier-Stokes solver based on immersed boundary methods [47,48]. It uses a Cartesian octree grid that can be dynamically refined to get a higher resolution of the fluid-air interface and the flow close to immersed objects. The software has been previously used to e.g. efficiently simulate robot application of sealing material and adhesives on automotive industry [49][50][51][52].
In the simulations a computational domain consisting of´10 10 1cubic cells is used, where each side of the cell is 1.0 mm. The dynamic refinements of the fluid-air interface cells is applied four times, making the smallest cell size 0.062 5 mm.    figure 1 at the prescribed speed of -10 mm s 1 in X-and Y-direction.

Rheology model
The rheology of the bioinks in IPS IBOFlow is modelled by a linear PTT-model [53]. This is a full stress tensor based model which captures the viscoelastic behaviour of the solder paste. The derived constitutive equation for the viscoelastic stress can be written as where λ is the relaxation time of the material, η is the viscosity contribution to the elastic stress, t ij is the viscoelastic stress tensor, S ij is the strain rate tensor and F is the PTT-model stress function. The stress function takes the nonlinear elongational behaviour of the viscoelastic stresses into account and is written as where ε is a parameter related to the elongational behaviour of the flow. In order for the model to work, the material related parameters (λ and η) need to be determined. Oscillation frequency measurements shown in section 2.1.1 give data of both storage modulus, G′, and loss modulus, G″, over a wide range of frequencies, ω. Both G′ and G″ can also be calculated as a function of the material parameters λ and η. The storage modulus is calculated as and the loss modulus is calculated as The complex modulus, G * , is calculated as are used for fitting of the parameters against the provided experimental data.

Surface tension force model
The surface tension force is approximated by a body force across the fluid-air interface and is modelled using the continuum surface force method by Brackbill et al [54]. It is formulated as where κ is the curvature of the liquid-gas interface, σ is the value of the surface tension between the two fluids and ∇α is the gradient of the volume fraction. The curvature is approximated as the divergence of the interface normal, wheren is the interface normal, which is in turn approximated by the normalised volume fraction gradient a a =   | | ( ) n . 8

Quasi-dynamic contact angle model
The quasi-dynamic contact angle model takes into account if the three phase contact line, the line where fluid-fluid-solid meet, is advancing or receding. Since there only is dependence on the contact line direction, two basic equations govern this model where  V cl is the contact line velocity, θ is the resulting dynamic contact angle from the model, θ A is the advancing contact angle and θ R is the receding contact angle. The contact angle is implemented as a boundary condition for the continuum surface force method by altering the interface curvature at the three phase contact line [55]. The interface normal is otherwise calculated according to equation (8), but to alter the interface curvature at the three phase contact line, it is instead calculated as cos , 10 f s t wheren s is the normal pointing outward from the solid surface andn t is the tangential surface normal pointing along the surface, towards the centre of the heavier fluid.

Rheology
Both of the tested bioinks, ink 6040 and 4% CNF ink showed a solid like behaviour in the oscillation frequency sweeps since the storage modulus was higher than the loss modulus (figure 2). Ink 6040 had a lower moduli and was therefore expected to flow more easily than 4% CNF ink.
In figure 3 the experimental data of G * and the resulting models of the parameter fitting are presented. The resulting values of λ and η for the respective bioinks are summarised in table 2.

Comparison of simulations and experiments
The rheology models of the complex modulus resembled the experimental data well enough to be used in the simulations. The tuned viscoelastic rheology model is then used in the framework for a 3D bioprinting simulation of the grid structure shown in figure 1. The same flow rate, printing speed, nozzle diameter and layer height is used in the simulation as for the experimental setup described in section 2.2. The results from the simulation are then compared to photos of experimentally printed grid structures. By using 4% CNF ink and ink 6040, which exhibit clear differences in viscoelastic properties, see figure 2, the ability of the framework to simulate different types of bioink was tested. The comparison for the 4% CNF ink, presented in figure 4, show close similarities between simulation and experiment. The uneven curvature of the printed line, as well as the thickening of the printed line before certain corners is accurately captured. The overhanging lines are quite distinct for the more viscous 4% CNF ink, which is seen in the simulations as well.
For the less viscous bioink, ink 6040, the flow is more driven by surface tension forces. Therefore the surface tension model and the static contact angle model plays an important role to accurately describe the interface curvature and surface flow. And as seen in figure 5, the smoother surface of the less viscous ink is reproduced in the simulations by the surface tension model. The static contact angle model correctly models the surface flow by governing the motion of the three phase contact line on the substrate. Overall the agreement between experiment and simulation is good, and this can especially be seen on the rounding of corners. Both the experiment and the simulation yields less distinct overhanging lines compared to the more viscous 4% CNF ink. Figure 6 shows the average height and width of the printed line compared to the simulated line. The simulation continuously measures the height and width of the first segment of the grid structure, while the experimental data is from multiple measurement points along the printed line. The standard deviation together with minimum and maximum values are also presented distinguish differences between simulation and experiment more clearly. The results from the two experiments verify the accuracy of the numerical simulation and indicate a well functioning rheology model. The simulation of 4% CNF ink has very accurate mean values for both line height and line width.    comes to printed line width. However, the overall comparison, including both the visual and the measurements, show that the simulation framework is able to capture the characteristics of ink 6040 and produce accurate results. Knowledge of the height and the width of a line printed with a certain bioink is important for determining suitable printing settings. The height will govern the suitable thickness of each 3D printed layer. If the thickness is set higher than the actual printed line height, the nozzle will after a few layers loose contact with the previously printed line causing the bioink to be dispensed inconsistently. In contrast, a too low thickness will cause layers being dispensed into each other, or alternatively hinder the flow of bioink out from the nozzle. The width is used when setting the distance between each printed line. If the intention is to print interconnecting lines to form a printed sheet, setting a distance larger than the line width will print lines with spacing in between while a too low distance will print lines that overlap each other. Thus, the simulation tells us that for printing 4% CNF ink with the flow parameters provided in section 2.2 (10 mm s −1 , 100 mg min −1 ), the layer thickness and line width should be set to m 400 m and m 420 m, respectively. On the other hand, for ink 6040 the layer thickness should be set to m 230 m and the line width should be set to m 780 m.

Varying printing parameters
The simulation showed good correlation with the experimental data and was therefore further used to  study how the printed lines of the bioink are affected by changing the printing speed and nozzle height.

Different printing speeds
The resulting layer thickness and line width presented in the previous section, are valid for a printing speed of -10 mm s 1 and a flow rate of 90-100 mg min −1 . Assuming a constant flow is kept, with increasing printing speed the line width should theoretically decrease, due to the lower volume dispensed per time unit. This correlation has also been shown experimentally in studies on 3D bioprinting where the line width controls the printing resolution [56]. By simulating, the change in line height, as a result of varying printing speed, can be predicted in addition to line width. Simulation can thereby be used to study how the cross section changes depending on the ink and the printing speed as shown in figures 7 and 8 for ink 6040 and 4% CNF ink.
The ratio of width to height, does not vary proportionally when changing the printing speed. For 4% CNF ink, the ratio calculated from the average line width and line height decreases from 1.5 to 0.9 when increasing the printing speed, and the data for both height and width is shown in figure 9. This result show the difficulty in experimentally determining the optimal layer thickness and line distance to set while printing. Also the standard deviation and min/max-values of the line width and the line height gives information on the printability of the inks at different printing speeds. For ink 6040, the standard deviation in figure 9 is larger at 5 and -20 mm s 1 indicating that the printed lines will have a more consistent shape throughout the printing at -10 mm s 1 and thus result in better printability. The prediction for the best suited printing speed of 4% CNF ink is more inconclusive, but it probably lies somewhere in between 10 and -20 mm s 1 . This is because the standard deviation of   The inks are intended for printing with cells for construction of synthetic tissues. There are several factors in the bioprinting process which may influence the cell survival; mixing cells with the bioink, osmolarity of the bioink, shear forces within the nozzle, diameter of the nozzle outlet, temperature, and, access to nutrients and oxygen. While dispensing the bioink, which is the simulated bioprinting step, there is a transition point from flowing vertically out from the nozzle to being dragged along the fixed printing plate. During this transition there is a build up of viscoelastic stresses on the outside of the nozzle tip, which may influence the cell survival. Therefore it is valuable that simulation can be used to predict these stresses, as shown in figure 10 where the distribution of viscoelastic stresses is visualised within a printed line. For both inks, the stresses are highest closest to the nozzle outlet and decrease further from the outlet along the printed path. Comparing the viscoelastic stresses between the two bioinks, the stresses are higher in the 4% CNF ink. This is an important factor for the characteristics of a 4% CNF ink printed line, whereas the surface tension forces are of greater importance for the characteristics of a ink 6040 printed line.

Different nozzle heights
The nozzle height is the distance between the outlet of the printer head nozzle and the printing plate upon which ink is being dispensed. The measurements of height and width of simulated printed lines is presented in figure 11. It is seen that decreasing or increasing the nozzle height with 0.1 mm, to 0.3 mm and 0.5 mm, respectively, does influence the line resolution, but not to a very great extent. For 4% CNF ink there is a trend towards better line resolution with  fewer fluctuations as the nozzle height is increased. Both the standard deviation and the difference between the minimum and maximum values are very low for a nozzle height of 0.5 mm. The same cannot be said for ink 6040, where the results suggest that using a lower nozzle height of either 0.3 or 0.4 mm would be favourable. At the highest nozzle height, the difference between the minimum and the maximum value becomes considerably larger, and an increase in standard deviation is observed as well.
Another factor which the nozzle height may influence is the printability of bioinks containing cells, since the distribution of viscoelastic stresses will change as shown by figure 12. For both inks, the viscoelastic stresses are clearly decreasing with increasing nozzle height, although a stronger trend is seen for 4% CNF ink. In the way of viewing changes to both nozzle height and printing speeds as tools to increase the printability of bioinks containing cells, by using them to reduce the viscoelastic stresses, the printing speed showed to have greater effect on the printability than the nozzle height.

Conclusion
As the field of 3D bioprinting advances, the need for simulations to predict the printability of bioinks and their influence on cells is increasing. Simulation will therefore be an essential tool for future development of 3D printed tissues. We have shown that CFD simulation with IPS IBOFlow is successful in capturing the dispensing of inks when 3D bioprinting and predicting the printing process for the studied parameters. The ability to simulate 3D bioprinting of inks with distinct difference in viscoelastic properties showed that the simulation tool can be utilised for a wide range of viscoelastic solutions and hydrogels. Furthermore, simulation is useful for determining the  impact on line resolution and viscoelastic stresses when changing printing parameters such as printing speed and nozzle height. Understanding which parameters affect the line resolution is vital for the mechanical properties of the printed structure. The ability to determine which printing parameters to change to reduce the effect of viscoelastic stresses on biochemical and biomechanical signalling is important for cell viability. In this study, simulation of a specific part of the 3D bioprinting process was presented: the flow of bioink out from the nozzle dispensed on the printing plate. Future development of the model used for simulating 3D bioprinting is expected to also consider the flow inside the nozzle to enable simulation of the ink from the point it starts to flow until the whole shape has been printed. Since simulation can aid in reducing costs, material usage, and time, it will be a valuable tool when developing and evaluating bioinks for 3D bioprinting in the future.

A.1. Lua script
The 3D bioprinting simulation in IPS IBOFlow is controlled through a Lua script. The first part of the script contains the setup of the simulation, which involves reading options from a separate options file, importing nozzle definitions and defining the boundary conditions used in the simulation. The boundary condition used for all boundaries, except the negative Z boundary, is an outlet boundary condition. In the negative Z direction a no-slip wall boundary condition is utilised, and both advancing and receding contact angles are prescribed at this boundary for the dynamic contact angle model. The path of the nozzle and the speed at which it follows it, is setup after the nozzle definitions has been imported, and it consists of multiple position points, which describes the printing path, and the nozzle speed. After these steps, the simulation is setup. The last part of the script contains the initialisation of the nozzle, activation of dynamic grid refinements and execution of the simulation.

A.2. Nozzle definition
The nozzle definition is kept in a separate file and read by the main script. It contains definitions of both the nozzle and the bioink. The nozzle definitions are shape, which is usually a cylinder, and diameter, which varies from nozzle to nozzle. The flow out of the nozzle is also defined, as either mass flow rate or volumetric flow rate. The bioink is defined by entering density, rheology model and parameters of the rheology model. A viscoelastic rheology model is used for 3D bioprinting and the parameters of the model are explained in section 2.3.1.

A.3. Options file
The options related to the simulation is stored in an options file. It contains options related to time settings, such as time step length, simulations end time and how often to write VTK output, which contains the data of the simulation at that time step. It also contains options for the surrounding fluid in the computational domain, which is the density and viscosity of air. VOF related options such as surface tension and type of dynamic contact angle model is also set here. Options related to the computational domain is also found here, by defining the number of base grid cells in each direction, and the length of a base grid cell.