Engineering considerations on extrusion-based bioprinting: interactions of material behavior, mechanical forces and cells in the printing needle

Systematic analysis of the extrusion process in 3D bioprinting is mandatory for process optimization concerning production speed, shape fidelity of the 3D construct and cell viability. In this study, we applied numerical and analytical modeling to describe the fluid flow inside the printing head based on a Herschel–Bulkley model. The presented analytical calculation method nicely reproduces the results of Computational Fluid Dynamics simulation concerning pressure drop over the printing head and maximal shear parameters at the outlet. An approach with dimensionless flow parameter enables the user to adapt rheological characteristics of a bioink, the printing pressure and needle diameter with regard to processing time, shear sensitivity of the integrated cells, shape fidelity and strand dimension. Bioinks consist of a blend of polymers and cells, which lead to a complex fluid behavior. In the present study, a bioink containing alginate, methylcellulose and agarose (AMA) was used as experimental model to compare the calculated with the experimental pressure gradient. With cultures of an immortalized human mesenchymal stem cell line and plant cells (basil) it was tested how cells influence the flow and how mechanical forces inside the printing needle affect cell viability. Influences on both sides increased with cell (aggregation) size as well as a less spherical shape. This study contributes to a systematic description of the extrusion-based bioprinting process and introduces a general strategy for process design, transferable to other bioinks.


Introduction
The technology of three-dimensional (3D) bioprinting provides the opportunity to create pre-designed, volumetric objects with a spatially defined distribution of embedded cells, enzymes and biofactors, which makes it a unique tool for diverse applications. The main focus in research has been so far on Tissue Engineering, but latest studies also have looked into applications in white biotechnology and the food sector; accordingly, mammalian cells, microorganisms and plant cells have been used for bioprinting [1][2][3][4].
To create macroscopic, volumetric structures in an appropriate time window and with high shape accuracy, extrusion-based bioprinting is the most suitable 3D printing technique [5,6]. A cytocompatible hydrogel of sufficiently high viscosity with integrated cells, a so-called bioink, is strandwise deposited and thus printed layer by layer. Challenges are the accomplishment of a high accuracy in relation to the computeraided design (CAD) model, a sufficient shape fidelity and stability over cultivation time, but also to retain high viability and specific activity of the cells [7]. This requires an ongoing process of material and process development. With regard to printing performance, fast printing with a high mass flow rate is desirable but leads to higher shear stress on integrated cells in the bioink [8]. Thus, an important aim is the identification of a critical threshold up to which the cells survive printing without severe damage. The effect of mechanical forces in printing on cell viability was found to depend on shear stress and exposure time determined by the geometry and length of the printing needle, the rheological behavior of the bioink and the flow rate [9,10]. The shear sensitivity of a cell increases with its size and decreases in case of a more spherical shape, while it is depending on the species and cell type, as well as on the constitution of the cytoskeleton and plasma membrane [11][12][13]. In this study, two cell types with a very different morphology have been compared: table 1 shows the main properties of a human mesenchymal stem cell (MSC) culture compared to a plant cell culture of basil. In contrast to MSC, the plant cell culture contains bigger, agglomerated cells with large vacuoles, rigid cell walls and higher variation in cell shapes. Therefore, the impact of shear stress is assumed to be higher for plant cell cultures than for MSC [14], but has not been studied in bioprinting so far.
A simulation of the material flow in bioprinting processes helps to understand interrelations between needle geometry, printing parameters, material properties and mechanical forces during extrusion. It can help to predict the mechanical stress on the cells or the shape of the extruded strand. Simulations can contribute to the optimization of the printing process and the establishment of new bioink compositions with less experimental work. A few studies have been published in the last years dealing with this topic regarding viability and physiology of different mammalian cell species in relation to the applied pressure and printing needle geometry [9,[15][16][17][18][19]. Increasing shear stress as well as residence time inside the needle was observed to affect cell viability [10,15,18,20]. In consequence, the flow rate is the most determining parameter for shear forces in the fluid flow and therefore the target value for cell viability considerations. In contrast, the indication of the printing pressure does not allow a comparison of the shear effects inside the printing needle for different bioinks.
Some recent publications address various methods for process characterization aiming at a better comparability and systematic process optimization [18,19,21,22]. For instance, Paxton et al present an approach that enables the general estimation of the 'printing window', the range of process parameters, for different bioinks which is limited by the printers pressure capabilities and the printing speed range feasible for printing (1-40 mm s −1 ) [18]. The present study focusses on a more general evaluation of relations between material properties, flow parameters and needle diameter which can serve as a valuable tool for systematic process optimization. Figure 1 shows the integrated approach of our study, which combines computational modeling as well as analytical calculation of the bioink flow with experimental validation and an investigation of factors influencing material flow, such as the cell density inside the bioink.
The bioink AMA containing alginate, methylcellulose and pre-gelled agarose in deionized water [3] was used as an example for pasty bioinks with complex rheological behavior with a view to meet the mentioned requirements in bioprinting.
Based on the Herschel-Bulkley law, describing non-Newtonian fluid behavior, a Computational Fluid Dynamics (CFD) simulation was conducted to analyze the material flow in the printing head consisting of a cartridge and a conical printing needle for different applied pressures and needle outlet diameters. Results of an analytical approach were compared with the simulation results in order to figure out if and to which extent the elaborate computational simulation is replaceable by an analytical calculation. The obtained flow profile was transferred to a dimensionless form enabling a fast characterization of the flow for different printing conditions. The results were summarized in a bioink-specific nomogram which allows the user to easily select printing parameters for an optimum set of values for printing speed, strand width, maximum shear stress and shear exposure time. The approach described herein for AMA can be easily applied to other bioinks.   of suspension carefully mixed into 1 g of the AMA bioink to a final concentration of 10 7 cells g −1 bioink.

Plant cell culture
An in vitro cell culture of basil (Ocimum basilicum L. var. purpurascens Benth., 'Cinnamon Basil') was used for plant cell bioprinting as described in [3]. The suspension culture was maintained in 50 ml Murashige and Skoog medium (MSmedium, including vitamins; Duchefa Biochemie B.V., Haarlem, Netherlands) with passaging every week by transferring 20% v/v of cell suspension to 80% v/v fresh MS medium. For rheological tests and bioprinting, plant cells were collected from MS medium at day 7 after passaging by using a glass frit filter (pore size 160-250 μm, ROBU ® GmbH, Hattert, Germany). Biomass was gently mixed into the prepared hydrogel bioink. The biomass dry weight content was estimated by gravimetric analysis of reference samples (n=10) after drying at 60°C for 48 h.

Preparation of the model bioinks
The bioink AMA containing alginate (30 g l −1 ), methylcellulose (30 g l −1 ) and agarose (10 g l −1 ) was used as described in [3] but it was prepared applying a slightly modified preparation method. Methylcellulose powder (Sigma Aldrich, Taufkirchen, Germany) was dissolved in deionized water at a temperature of 4°C. Alginic acid sodium salt (order number: 71238, Sigma Aldrich, Taufkirchen, Germany) and agarose (UltraPure™, Thermo Fisher Scientific, Waltham, USA) were added, dissolved while heating to 90°C and autoclaved at 121°C for 20 min. While cooling down to room temperature the AMA-blend was stirred to obtain a homogenous bioink and allowed to rest at 4°C at least overnight. The second blend for experimental validation of flow simulation contained only 30 g l −1 methylcellulose (M-30) prepared equally to AMA.

Characterization of the model bioinks
Using a syringe, a defined volume of cell-free bioink was weighed to calculate the density r (n=10). Rheological properties of AMA and M-30 were evaluated at 20°C using a rotary rheometer (Rheo-testRN 4.1, Medingen, Germany) with a plate/plate measuring system (Ø30 mm) and a gap of 0.1 mm.
Shear thinning experiments were performed by applying a linearly increasing shear rate from 1 to 600 s −1 (increments of 0.08 s −1 ) to obtain the corresponding shear stress. To determine the yield stress t 0 of AMA, a stress-controlled measurement of the shear rate was performed in the range of  Pa. An almost constant shear rate indicates an elastic deformation of the tested material. A sudden increase of the shear rate is related to a break down (yield) of the elastic structure. The applied shear stress at this point is defined as the yield stress t . 0 This method was selected out of many existing approaches to determine the yield stress as described in [24]. For normal force measuring to determine first normal-stress difference N , 1 the rotary rheometer MCR 300 (Anton Paar, Germany) was used with an adapted gap of 0.5 mm in a range of the strain rate from 1 to 200 s −1 . All experiments were conducted three times.

3D bioprinting (extrusion)
The Bioscaffolder 3.1 (GeSiM mbH, Radeberg, Germany) was used for extrusion-based 3D printing. The AMA bioink was transferred into a 10 ml cartridge (Nordson Corporation, Westlake, USA) and dispensed through conic dosing needles (Globa-coGmbH, Roedermark, Germany). To validate simulation, the mass flow rate m  was determined by gravimetrical analysis of the amount of dispensed bioink for different applied pressures (from 40 to 300 kPa) and needle outlet diameters (0.25, 0.41, 0.61, 0.84 and 1.2 mm), each for a time of 20 s. The distance between the needle tip and the petri dish was about 100 mm to avoid accumulation of the bioink at the outlet. In cell stress tests, printing pressures of 150 and 300 kPa and outlet diameters of 0.25 and 0.84 mm were applied to compare moderate and extremal shear conditions. All experiments were conducted three times.

Cell viability test
Cell viability was determined via cell staining in the suspension, in the bioink before and directly after printing as well as after rheometer-mediated shear treatment. All dyes were obtained from Invitrogen by Thermo Fisher Scientific, Carlsbad, USA.  Figure 2(a) shows the geometrical dimensions of the 2D model. The geometry of the printing head was taken from two 3D-CAD files kindly provided by the manufacturers (cartridge: Nordson Corp. Westlake, USA; needle: Globaco GmbH, Rödermark, Germany). The assembly was merged and simplified by replacing the setback between needle and the cartridge with a 45°edge in order to generate a homogenous mesh and to obtain a stable solution. The needle outlet diameter D was changed by slicing the geometry with a line through the fix point A and removing excess material. The cartridge had in real an overall length of 10 cm, the inlet section of the domain was fixed to a length of approx. 15 mm to ensure smooth inlet flow conditions and no influences on the distribution of flow parameters at the outlet. The outlet had no further extension. The structured mesh (figure 2(b)) was generated with Ansys Meshing 18.2 and consists of 3600 quadrilateral elements (quads) with middle node and a maximum element size of 0.5 mm. The cell size of the 30 elements in radial direction was refined towards the wall with a smooth geometrical transition factor (spacing) of 1.1.
Ansys Fluent 18.2 was used to perform a steadystate 2D-axisymmetric calculation of the printing process. The flow was modeled isothermal and laminar, assuming constant density (see 4.1), using a Herschel-Bulkley fluid model according to where t is the shear stress, t 0 is the yield stress, g  the shear rate, n the flow index and k the consistency factor.
The boundary conditions for the model are shown in figure 2 Dp meas was applied at the inlet, where the ambient (reference) pressure p 0 of the experiment was assumed constant at 100 kPa while Dp meas was the measured overpressure in the experimental setup. Rather than the ambient pressure p , 0 the measured mass flow rate m  obtained from the experiment was set as the outlet boundary condition since it is assumed to be the leading influence factor on the shear stress magnitude at this cross-section. Hence, the (area-averaged) static pressure p out at the outlet resulting from the simulation may differ from the software-adjusted pressure of 100 kPa (see below). While the centreline of the cartridge is the symmetry axis of the domain, all radial borders were modeled as stationary noslip walls with zero wall velocity ( = u 0  ). Calculation was performed using pressure-velocity coupling (SIMPLE or SIMPLE-C algorithm) and second-order discretization for pressure and momentum and least-squares cell-based discretization for gradients. Relaxation factors were 0.5 for pressure, 0.8 for momentum and 1.0 for density and body forces. With these settings, relative residuals of 10 -5 were reached after approximately 500-1000 iterations for all calculations which corresponds to a calculation time of only several minutes on a local machine with three processor cores. Area-weighted average values for the static pressure and the shear stress at the outlet were monitored approaching a state of convergence even faster, confirming that a stationary solution was obtained.
For further analysis regarding impact on the cells, the radial distribution of the velocity u, the shear rate g,  the dynamic viscosity h and the shear stress t xy at the outlet line as well as the axial discharge of the pressure p along the centreline (symmetry axis) were evaluated for the converged solution.

Analytical calculation
An analytical algorithm [25] based on the analytical solution of the flow of the Herschel-Bulkley fluid in arbitrarily shaped axisymmetric channels was used to predict the distribution of flow parameters in the printing needle even faster. The equations given in [25] are based on the assumption that the inclination angle of the walls are small in axial directions and, thus, the flow can be considered purely axial. Further, it can be proved that the given equations are also valid to be used with dimensional values. In contrast to the original paper, the geometry of the needle wall is not described by analytical functions but by discrete values of the radius R z ( ) over the axial coordinate z ( figure 3).
The equations in [25] have been solved at the midside nodes of the outer elements of the contour line: In consequence, the axial length of a contour element is then: The value of the axial pressure gradient =q z dp dz ( ) / was determined iteratively for each element midside radius R m i , using Newton-Raphson method and applying first-order upwind difference scheme for ¢ f q :  is the volume flow of the incompressible fluid through the cross section pR .
A mesh sensitivity study was performed to evaluate the influence of the number of axial divisions on the overall pressure difference. Therefore, the elements of the basic case (figure 3) have been further subdivided into = A 2, 3 or 4 additional elements or have been merged pairwise omitting every other node ( = A 0.5), (table 2). As can be seen in figure 4, the differences in the printing pressure are negligibly small and do not change significantly for  A 3. Hence, the 'mesh' with 367 axial divisions was used for further investigations.
After calculating the axial pressure gradient q out exactly for the outlet radius R out according to equation (5), the radial distribution of the flow parameters at the outlet can be determined analytically, i.e. no discretization in radial direction is necessary.
The radius and the velocity of the plug flow region at the outlet are calculated according to The radial velocity profile can then be described by which yields the following shear rate profile: The shear stress distribution follows from the Herschel-Bulkley law (1): At last, the dynamic viscosity over the radius can be calculated from The correlation coefficient for the curve fit is = R 0.98. 2 The low value of the flow index n of 0.2 indicates a distinct shear thinning behavior. The yield stress t 0 of the bioink was determined by a sudden increase of the shear rate, which indicates the limit of the linear elastic region (figure 6(b)). ( ) describing the axial contour of the printing needle. The given example corresponds to the outer nodes of the CFD model (see figure 2) and it is the reference case for the mesh sensitivity study, = A 1.

Validation of simulation and analytical calculation with experimental data
In order to validate the simulation and the calculation model, the overall pressure difference between the inlet and the outlet Dp was compared with the set values in the experiment for the different mass flow rates m  (figure 7, supplementary data is available online at stacks.iop.org/BF/12/025022/ mmedia). As figure 7(a) shows, the simulated and the calculated pressure differences were nearly identical. Both approximation curves can be described with a power function, where the exponent is close to the flow index of the used Herschel-Bulkley material model, = n 0.2. However, there is a significant systematic deviation from the experimentally measured curve (figures 7(a) and (b)). It is assumed that this can be explained mainly with elastic effects (section 4.4).     Figure 8 shows the distribution of velocity u, shear rate g,  shear stress t and viscosity h over the radius of the needle outlet obtained from CFD simulation and from analytical calculation. The flow shows a shear thinning behavior, which leads to a plug-like flow. In contrast to a fully developed, parabolic flow profile of a Newtonian fluid, the velocity in the needle center stays almost constant in the center of the needle area (plugflow region) and decreases fast near the wall [26].
The example in figure 8 demonstrates that the analytic approach serves well to predict the flow in the printing needle in comparison with the simulation. While differences in the axial pressure discharge ( figure 9) are negligibly small, the calculated velocity profile is flatter in the middle of the channel than the CFD profile ( figure 8(a)). According to the simplified definition in equation (12), this leads to a calculated shear rate of zero for = r 0, whereas the shear rate used in Ansys Fluent is defined to be the second invariant of the shear rate (rate-of-deformation) tensor = D D ij ( ) and, thus, contains also components of the velocity gradient in other directions:  (15) can be simplified to formula (12). However, since the printing needle has a small taper angle, the shear rate calculated by the CFD program has an offset compared to the analytical solution that is based on the assumption (12), figure 8(b). Consequently, the shear stress in the channel core yields higher values than the analytical approach, the latter of which gives the yield stress t = 60 Pa 0 at the axis. This results in finite dynamic viscosities in the center as opposed to equation (14), which calculates h  ¥ for  r 0.

Normalization to dimensionless flow parameters
As the dimensionless equations in [25] indicate, the problem allows a self-similar solution, referring all involved quantities to specific characteristic values. In the present case, we referred the velocity to the maximum velocity at the outlet, which is located at the centreline are also coincident ( figure 10(b)). With the Herschel-Bulkley law (13), we could also get a dimensionless that shows self-similarity ( figure 10(c)). To find a similar approach for the dynamic viscosity h t g = , /  we needed to make a simplification that the occurring shear stresses are much greater than the yield stress, t t, 0  which is true for all investigated CFD cases. When plotting over r , * the maximum deviation of the normalized viscosity h* is 3.9% and the standard deviation is only 0.14% for all cases. Hence, nearly no differences would be seen in figure 10(d), which displays the arithmetic average over all cases. For the analytical solution, there actually exists no valid curve to be plotted in figure 10(d), since h = ¥ max (see equation (14)). For comparison, the value of h max for each CFD-case (p , In m  ) was also taken as reference value for the analytical solution. We could also find self-similarity in the pressure discharge along the axis of the printing needle for constant needle diameter D (figure 11), if we normalize the pressure with the area-averaged values of the pressure at the inlet (p in ) and the outlet (p out ): and print it over the normalized z-axis where z In is the axial coordinate at the inlet and L is the distance between the inlet and the outlet. Figure 12 shows contour plots over the whole CFD domain for all reduced parameter described above.
All parameters are ranging from 0 to 1, except for the viscosity in figure 12(e), which also takes values h 1, *  especially in the regions with low shear, e.g. in the inlet section.   respectively. Despite of major deviations between the calculated and the simulated maximum shear rate ( figure 13(b)), the calculated values of the maximum shear stress t max (figure 13(c)) lay only 4.5%-5.7% above the simulated points and for the dynamic viscosity h max (figure 13(d)) the differences are even smaller.
Hence, the analytic approach is suitable to quickly predict the maximum shear stresses on the fluid in the outlet with reasonable precision. Hence, CFD simulation is not necessarily needed. However, the calculation is based on simplified assumptions regarding the mechanical stress in the center of the printing needle, so that simulation is still complimentary to map the mechanical forces over the complete printing head.
For the application in bioprinting, the effective printing speed (feed rate of the printing head, i.e. mean flow velocity at the needle outlet) lies typically between 5 and 25 mm s −1 which corresponds to a mass flow rate between 1 and 5 mg s −1 for a 0.61 mm needle. Shear thinning of the bioink leads to a smaller increase of shear stress with increasing mass flow rate in comparison to a Newtonian fluid. For example, the maximum shear stress is 2300 Pa for a mass flow rate of 2 mg s −1 and increases to 2632 Pa for 4 mg s −1 . While the maximum velocity increases linear to mass flow rate by 100%, the increase in shear stress is only 14%.

Process design tool for bioprinting
Calculation results were further used to create the nomogram in figure 14, which allows to quickly estimate the overall mass flow rate and the maximum shear stress at the outlet of a given printing needle with diameter D for an applied pressure Dp or any other combination of three parameters, e.g.
The self-similarity of the flow behavior allows to easily receive results for different mass flow rates and needle diameters just with the input of the flow curve and the density of the bioink of interest. A nomogram as in figure 14 can be created for different bioinks, which enables the user to find a bioink with suitable rheological characteristics as well as an optimum applied printing pressure and needle diameter with regard to processing time, shear sensitivity of the integrated cells, shape fidelity and strand dimension. The example in figure 14 shows that the maximum shear stress in a 0.25 mm needle (A) is only slightly lower than in a 0.84 mm needle (B) for the same applied pressure, while the mass flow rate in A is 97% lower than in B. Depending on the aspects of designing the printing process, the nomogram can also be adapted, i.e. to plot the residence time in the needle as a further influence factor on the cell viability or the mean flow velocity in order to assess the production speed.
The optimal printing velocity can be determined with the assumption that the printing speed should be equal to the mean extrusion velocity u. With    15(b)), confirming its mainly viscous behavior.
In figure 16, the experimental mass flow rate is plotted for different needle sizes and a constant applied pressure of 200 kPa for AMA or 40 kPa for M-30, respectively. As the mass flow rate of AMA increased with the printing needle diameter, which can be approximated by a power law, also the difference between calculated and experimental pressure increased. In contrast, the deviations between the pressure values from the experiment and the calculation were negligibly small for M-30, also at higher needle diameters, for which the mass flow rates increased linearly. As a conclusion, it is expected that also for other bioinks without pre-crosslinked components like gelatine methacryloyl (GelMa), elastic effects are negligible within the usual range of printing conditions and the results from simulation and calculation can be used without further correction. Elastic effects then appear only at higher flow rates or for extremely small needle diameters, both of which may not relevant for bioprinting.  In contrast, AMA contained alginate, highly disperse methylcellulose [27] and gelled agarose. Especially in converging flows, long non-crosslinked molecules are able to align in the direction of flow leading to significant differences between axial and radial normal stresses and, thus, to a pronounced anisotropic, viscoelastic behavior, which cannot be described solely by viscous laws as in the shown simulation or in the calculation [28].
Oldroyd et al assumed that in a simple steady-state shear flow, elastic effects mainly appear as first normal stress difference where s xx acts normal to the flow direction and s yy in flow direction. The ratio between elastic and viscous forces in a viscoelastic fluid is characterized by the Weissenberg number Wi [28], which can be expressed according to the upper-convected Maxwell model presented by Oldroyd [29,30] as: By measuring g N 1 ( )  and t g xy ( )  simultaneously (figure 7(a)), it is not necessary to determine a relaxation time l of the bioink in order to get Wi.
The calculated shearing profile directly at the needle outlet is still expected to correspond with the real shearing conditions. Reasonable since normal stress effects are assumed to affect mainly the flow parameters in axial direction and the measured flow rate was set as outlet boundary condition. The local Weissenberg number according to equation (23) was exemplary calculated for the maximum shear rate g R ( )  at the outlet. Wi increased linear with increasing printing pressure ( figure 15(a)), while also the relative deviation between experimental and simulated pressure showed a similar tendency ( figure 7(b)). The pressure loss due to the elastic effects in the fluid flow was therefore assumed to be the main reason for the deviation of the calculation results from the experiment. In a following study, a correction should be integrated in an extended model to predict the mass flow rate for an applied pressure.
It has to be proven if the simple approach only considering first normal stress difference is valid for highly elastic bioinks or if more complex software, like Ansys Polyflow ® , 3D simulation and multi-mode models are necessary [31]. Complex fluid models may serve to model flow behavior of viscoelastic fluids in a more realistic way for Weissenberg numbers > Wi 1 and a wider relaxation time spectrum like in [31,32]. Anyway, the rheological measurement of g N 1 ( )  besides t g xy ( )  gives a first idea of the relevance of elastic effects for printing with the respective bioink. Finally, the model was validated exemplarily for different needle outlet diameters at a constant printing pressure of D = p 200 kPa.

Rheological effects on pressure drop
The validity of both the simulation and calculation model is limited by the underlying rheological data. For small needles and high mass flow rates, the shear rate exceeds the maximum value of 600 s −1 that can be measured with the used rheometer. Higher shear rates led to an ejection of the bioink out of the rheometer gap, measuring even higher shear rates would require a capillary rheometer. The model in Leppinemi et al also overestimates the mass flow rate with a pasty bioink containing nanofibrillated cellulose and alginate, which yields a similar flow index as AMA. They explained this behavior with a plateau of shear stress at shear rates extending the measured range, which leads to a deviation from the used power law model [33]. This effect could also appear in our experiments. Anyway, we already observed a deviation of the model from the experiment at mass flow rates that resulted in shear rates, which were definitely in the measured range, confirming the influences of elastic effects. Other studies just calculate the flow inside the printing needle without considering the cartridge [15,18,19]. As shown in figure 12(a) for our system, in the narrowing part of the cartridge appears already a pressure drop of about 20% of the overall pressure head.
The calculation in [18] also lead to an overestimation of the mass flow but the authors of this work suspected wall-slip as possible reason. When varying the gap of the rheometer between 0.1 and 1 mm for AMA, no change in the measured flow curves was observed indicating the absence of wall-slip (data not shown). However, it is assumed that the experimental flow rate with the presence of wall slip would be even higher than the calculated value rather than lower.

Cell concentration and aggregation
For calculation of the fluid flow in the printing process, it is important to study the impact of cell concentration, size, shape and aggregation on the macrorheological behavior of the bioink. In figure 17, the flow curves of AMA with different concentrations of hTERT-MSC (a) and different concentrations and aggregation sizes of basil cell culture (b)-(c) are shown.
Human cell culture. In the presence of hTERT-MSC, the viscosity did not change significantly. For a high cell concentration of 20×10 6 cells g −1 bioink, the difference was marginal compared to the bioink without cells. As expected, no significant impact on the rheology was observed for cell densities usually used for bioprinting, because the volume fraction for the highest cell number of 20×10 6 cells ml −1 with a cell diameter of 40 μm was still less than 150 ppm. Furthermore, cells have a water content of >90% and, thus, nearly the same density as the bioink, and they are to a certain extent elastically deformable [11]. Due to the small cell diameter of less than 10% of the needle outlet, a narrow particle size distribution and the spherical shape of the hTERT-MSC, the macroscopic fluid behavior can be considered continuous [34], which confirms the application of the shown simulation model.
Also other studies revealed no impact on the rheology of the respective bioink for similar cell concentrations: hepatocytes in a similar concentration as in the present study of 10 7 cells ml −1 bioink showed no influence on the rheological behavior [35]. Manojlovic et al tested even higher cell concentrations of 5×10 8 and 1×10 9 cells ml −1 with yeast, which did not affect the rheological properties of non-crosslinked 20 g l −1 alginate [36] either. A slight increase of viscosity and decrease of the flow index was found for a less viscous bioink of 15 g l −1 alginate when fibroblasts were added to a final concentration of 1 and 10×10 6 cells ml −1 [37]. It should be noted that the cell-laden bioink was only compared with the initial, non-diluted bioink; a reference without cells but with an equal amount of culture medium was not reported [37]. In contrast to our approach, that approach does not allow to observe the effect of cells separately from dilution and change in ion concentration when adding culture medium to the bioink.
Plant cell culture. For the bigger, aggregated plant cells enclosed by a rigid cell wall, the cell concentration in the bioink had a considerable effect on the rheological behavior of the bioink. A pronounced peak in the viscosity curve below a shear rate of 100 s −1 was observed which increases with the cell agglomeration size ( figure 17(c)). As confirmed by microscopy images (see 4.5), the study in the rheometer lead to a destruction of bigger aggregates at the beginning of the shearing due to the small gap between the plates of 100 μm, which interfere with the measurement. The cell diameter should be significantly smaller than the needle diameter. Otherwise, a continuous simulation is not valid and the cells have to be considered as a second phase.
At higher shear rates between 100 and 500 s −1 , the viscosity of AMA containing different cell concentration and aggregation size converged to the same values. At 500 s −1 , a slightly higher viscosity of about 1.5 Pa s was observed for cell-laden samples, which was not proven systematic. For plant cell suspensions in culture medium, it is known that the viscosity is strongly dependent on the biomass concentration [38,39]. The variety of cell morphology, aggregation and cell size depending mainly on the water content of the cells interfere with the effect of biomass concentration [12]. For simulation, the flow curve needs to be measured for various cell concentrations and with a gap size in the same range as the diameter of the used printing needle.

Cell viability after shearing/printing
Cell viability of hTERT-MSC and basil cell culture embedded in AMA was observed directly after printing as well as after treatment in the rheometer in order to study the influence of the different shearing mechanisms on the cells. The results are shown with the corresponding maximum shear stress and residence time out of the analytical calculation results.
A wide range of flow conditions was tested to assess the impact of shearing on the cell viability during the printing process, including applied pressure of 150 and 300 kPa and needles with outlet diameters of 0.25 and 0.84 mm. Alternatively, a constant shear rate of 10 s −1 , which correspond to an average shear stress of t = 1312 Pa xy (figure 6(a)), and 800 s −1 (t = 3068 Pa xy ) was applied with the rheometer, for an exposure time of 10 s (similar to the printing process) and significantly longer time of 300 s. The ratio of living cells in the initial cell suspension and in the bioink before processing serve as reference values.
Human cell culture. In figure 18, the live cell ratio of hTERT-MSC is shown. Viability after printing with Figure 18. Viability of hTERT-MSC in AMA after extrusion printing and shearing in the rheometer. The live cell ratio was determined for the initial cell suspension in medium ('suspension'), cells mixed into AMA before ('bioink') and directly after printing ('printed bioink') as well as directly after shearing in rheometer. For printing, the applied pressure Dp and needle diameter D were varied. For the measured mass flow rate m meas  the mean residence time t res and maximum shear stress at the needle outlet t max was calculated. The shear rate g  in the rheometer varied with the calculated uniform shear stress t and duration t. (mean±SD of n=9, significance ( * p>0,05).) a 0.25mm needle was on average 48% and, thus, slightly lower than with the 0.84mm needle (about 60%). Although the residence time t res and calculated shear stress t xy, max were considerably different the tested printing conditions did not lead to significant differences in live cell ratio directly after printing compared to the non-printed reference bioink. Even with the 0.25 mm needle and a mass flow rate of 190 mg s −1 , which is significantly higher than for usual printing conditions, the decrease of living cells was not significant. Only in rheometer test, the live cell ratio was reduced.
Several authors reported a decrease of the number of living cells after printing by 20%-50%-the smaller the needle diameter and the higher the applied pressure was, the higher was the loss in cell viability [8,9,40]. Contrary to our study, all of them used cylindrical printing needles and less viscous alginate solutions. Li et al compared cylindrical with conical needles regarding the mass flow rate and observed a reduction of live cell ratio of less than 5% for conical needles and resulting in a drop of up to 20% for cylindrical needles [8]. A cylindrical needle causes a long-lasting high shear stress contrary to conical needles, where the exposure to extreme shear stress is limited to the lowest part of the needle near the outlet, as shown in figure 12(d).
In rheometer treatment, the shear stress (3068 Pa) was chosen as similar to the maximum shear stress when printing with a 0.84mm needle and 150 kPa working pressure. Therefore, the significantly lower live cell ratio of the rheometer-sheared hTERT-MSC can be explained by the longer exposure time. While in rheometer the cells are permanently exposed to the indicated shear stress, shearing varies in printing for each cell depending on the pathline of the particle. t res indicates the average time one particle needs to pass the part of the printing head (35 mm), where shear stress mainly appears. Other authors identified a relation between residence time and viability earlier [15,18]. According to Li et al an increasing length of the cylindrical needle leads to a larger decrease of surviving cells [15].
The flow index of AMA of n=0.2 in the present study was significantly lower than for the bioinks in the mentioned studies where n>0.5 indicating a stronger shear thinning behavior of AMA. A lower flow index leads to a more plug-like flow profile [26]. For AMA, the simulated shear rate was in 55% of the needle outlet area smaller than 20% of the maximum shear rate at the wall ( figure 10(b)). Li et al [22] expected that a plug-like flow of a highly shear thinning bioink leads to a protection of the cells against shearing. Anyway, comparison between different studies is limited by different bioink rheology and sometimes missing information about flow rate data, printing head geometry and simulation data. It is important to mention that the flow rate as directly shear-determining parameter is more reliable to compare shear conditions for needle diameters and bioinks than the applied pressure.
Cells behave viscoelastic. To damage a cell mechanically by short-term load, a combination of large tangential deformation and high shear rates are needed. Small and sudden deformation leads to an elastic response of the cell while slow and high deformations can be compensated to a certain extent by cell elongation [11]. Critical shear stress still needs to be tested with regard to cell type, shape and physiology and it is greatly influenced by lipid migration within the cell membrane caused by mechanical deformation [11]. Anyway, the small cell diameter of less than 10% of the needle outlet and the spherical shape in the bioink of the tested MSC minimize the effective shear stress. In extrusion printing with hTERT-MSC and AMA the critical shear stress leading to cell dead is not reached for a wide range of needle diameters and flow rates. Therefore, this aspect does not need to be considered in process optimization.
Compared to the initial cell suspension with a live cell ratio of 90%, the viability was significantly lower for the cell-laden bioink. Mixing the cell suspension into the bioink was crucial for a significant decrease in viability for hTERT-MSC. So it seems that other stress factors in fabrication process were more critical than the shearing inside the printing needle including the mechanical stress due to mixing, a dense polymer network, static pressure and maybe osmotic stress in the bioink.
Plant cell culture. In figure 19, images of live-deadstaining of basil cell culture are shown. Due to aggregation and heterogeneous shape and size of the cells, a quantitative analysis was not possible. However, viability was compared between the cell suspension and the cell-laden bioink prior to printing ( figure 19(a)). After printing, a higher amount of dead cells was observed in the printed bioink compared to the non-printed bioink, although the majority of the cells were still alive ( figure 19(b)). The higher impact on cell viability compared to hTERT-MSC is not surprising as for cell cultures with bigger cells of more various shape, effective mechanical forces on a single cell are higher than for small spherical cells exposed to the same shear stress [41,42].
The application of bioprinting to bioprocesses with plant cell cultures is a new approach introduced by us recently [3]. To date, no other study investigating mechanical stress on plant cells in bioprinting has been reported. In vitro cultivation of plant cells has great potential for sustainable, bio-based manufacturing of active substances for pharmaceuticals, food and cosmetics, as GMP-compliant or continuous [43]. Bioprinting allows a controlled immobilization in matrices, free-formed and adapted to the bioreactor for optimal mass transfer [3]. The cells are protected against shearing during cultivation but exposed to a short, intense shear force in printing. Many studies reveal effects of fluid forces on plant cell cultures, which are commonly cultivated in suspension, depending on the shear stress and exposure time. They aimed to study long-term effects of stirring or aeration in cultivations [42,44,45]. It was found, that nonacute loading affect the cells over the whole cultivation [44] and shear stresses already 20 times lower than for bioprinting with AMA are critical for cell viability. This is supported by our rheometer experiments, which revealed a decrease in viability in samples treated at a shear rate of 10 s −1 only after 300 s and of 800 s −1 already after 10 s ( figure 19(c)). A shear rate of 10 s −1 for 10 s did not led to a visible effect on the viability. The corresponding shear stress is 1300 Pa, which lead to following conservative assumption for the printing process: The maximum shear stress in printing should not exceed 1300 Pa. Regarding the nomogram in figure 14 this means for printing with AMA and a 0.84 mm needle that the applied pressure should not exceed 50 kPa which results in a mass flow rate of 1.5 mg s −1 corresponding to a printing velocity of 2.75 mm s −1 .
At higher shear stress the exposure time plays a secondary role while cell destruction already appears earlier. Even though in the recent experiments there were still about 60%-80% living cells directly after printing, apoptosis and sub-lethal effects may appear [12] as reaction to mechanical stress. In a previous study, we showed that the majority of the cells survive printing and proliferating in AMA for several weeks [3].
Cell aggregates were disrupted already by mixing into the bioink ( figure 19(a)) and to a smaller extent later in the printing process ( figure 19(b)). The longer duration of shearing and the small gap in the rheometer led to even Figure 19. Viability of basil cell culture in AMA after extrusion printing and shearing in the rheometer. The live cell ratio was determined (a) for the initial cell suspension in medium ('suspension'), cells mixed into AMA before ('bioink'), (b) directly after printing ('printed bioink') and (c) directly after shearing in rheometer. For printing, the applied pressure and needle diameter were varied. For the measured mass flow rate m meas  the mean residence time t res and maximum shear stress at the needle outlet t xy, max was calculated. The shear rate in the rheometer varied with the calculated uniform shear stress t xy and durationt. smaller aggregates ( figure 19(c)). As for several plant species, aggregation of the cell cultures is needed to maintain metabolite production and growth [46], the effect of aggregate size on their performance has to be investigated in further studies.

Conclusion
The present study introduces an approach to describe the fluid flow of a bioink in the printing head during extrusion. Numerical simulation and analytical calculation were used to predict the mechanical stress on the cells, the pressure gradient and flow rate and hence, contribute to the optimization of the printing process and the establishment of new bioink compositions with less experimental work.
The analytical calculation method reproduced the results of the CFD simulation well concerning pressure drop over the cartridge and maximal shear parameters at the outlet. Introducing dimensionless flow parameters allowed finding a self-similar solution to describe the fluid flow, valid for different needle diameters and applied pressures.
On this basis, a bioink-specific nomogram was created, which enables the fast and easy identification of suitable printing parameters under consideration of the shear sensitivity of the cells, the strand diameter and the overall process performance. The detailed description of the calculation in the present study allows creating nomograms for different bioinks, without the need of a numerical simulation. This study contributes to a systematic description of the extrusion-based bioprinting process and should be combined for optimization with simulation of the strand deposition like in [21].
In experimental validation, the viscous flow model was tested to be sufficient for the description of the shearing conditions inside the printing head. For a bioink containing methylcellulose (M-30), which serves as example for a bioink without pre-gelled components, the calculated data fits well to the experimental data for a broad range of the flow rate. For the model bioink AMA the calculation overestimates the mass flow rate with increasing applied pressure. Rheometrical determination of the Weissenberg number Wi as ratio of elastic effects (first normal force difference (N 1 )) to viscous effects (shear stress t xy ) revealed that the elasticity of the bioink hampered the fluid flow already inside the printing device. Rheometrical measurement of N 1 besides shear stress is a simple general test that allows evaluating the dominance of elastic effects in the fluid flow.
Bioprinting with the bioink AMA did not affect the viability of the hTERT-MSC due to shear stress. Basil cell culture showed smaller aggregates and a lower viability after printing, which points out the necessity to better adapt the printing parameters on bigger (plant) cells and cell aggregates concerning flow rate and needle diameter. In a rheometer with long exposure to high shear rates, a shear stress and time-dependent decrease of the viability was observed, which revealed the exposure time as a crucial factor. Complex interrelations make a basic investigation on the biocompatibility of bioink and every step in fabrication for every cell type necessary.
Future work will focus on the characterization of the printing process for several bioinks with different cell types. Results of the flow calculation can be summarized in nomograms. Using the limiting parameter out of the cell shear tests in the rheometer, like maximum shear stress or residence time, optimal parameter, regarding process performance and cell viability, can be found. This systematic description is an effective tool to optimize the bioprinting process.