Bioink formulation
3D bioprinting bioinks were formulated using bovine gelatin (Sigma, France), low viscosity alginate (Sigma, France) or very low viscosity alginate (Alfa Aesar, Thermo Fisher France) and fibrinogen from bovine plasma (Sigma, France). Stock solution of 0.2 g/mL gelatin, 0.04 g/mL alginate and 0.08 g/mL fibrinogen were dissolved overnight at 37°C in DMEM (without calcium, with glutamax1, Invitrogen, France) supplemented with 10% foetal calf serum (HyClone, USA), 20 µg/ml gentamicin (Pantapharm, France), 100 UI/ml penicillin/streptomycin (Sarbach, France) and 1 µg/ml amphotericin B (Bristol Myers Squibb, France). Table 1 presents the composition of the different formulations used during this study.
Bioink rheological characterisation
A stresscontrolled rotational rheometer (DHR2, TA instruments, New Castle, USA) was used to evaluate the rheological behaviour, the viscoelastic balance and the static yield stress value at 21°C with a 25 mm plate geometry. For rheological behaviour (storage modulus G’), a shear rate sweep was applied to hydrogels from 0.01 to 10 s− 1. The viscoelastic balance was evaluated thanks to the phase angle by shear dynamic analysis at frequencies 0.015 Hz, considered to be at rest, and 1.5 Hz, considered to be at flow. Then, the yield stress value was defined by 3iTT test with a shear rate at 30 s− 1 (mimicking extrusion bioprinting) during 120 s. With this test, the yield stress was calculated with the SAOSLVER method.
3D Bioprinting procedure
Just before printing, each cell type was trypsinized, suspended in fibrinogen to obtain the chosen cells/bioink mL seeding (see Supplementary Table 1) and then formulated with gelatin and alginate to produce the final bioink. After homogenization, a 3 or 10 mL sterile cartridge (Nordson EFD, France) was filled with the bioink, incubated 15 minutes at 37°C and then 30 minutes at room temperature to stabilise the bioink rheological properties. The cartridge was then loaded in a 6axis robotic bioprinter (BioAssemblyBot®, Advanced Lifescience Solutions, USA) or a BIO X (CELLINK, Sweden) and used to print standardised 0.3 cm3 bioprinted tissues (1 cm*1 cm*3 mm). A 410 µm diameter, 6.35 mm long needle (Nordson EFD, France) was used to bioprint at a set speed of 10 mm/sec.
Table 1
Final composition of the 10 bioink formulations used in the present study.
CODE

Fibrinogen final g/mL

Alginate very low viscosity final g/mL

Gelatin final g/mL

Gelatin proportion g/total g

BI1

0.02

0.02

0.05

1.25

BI2

0.02

/

0.1

3.3

BI3

0.02

0.015

0.075

2

BI4

0.02

0.005

0.125

5

BI5

0.02

0.025

0.025

0.56

BI6

0.02

0.03

0

0

BI7

0.02

0

0.15

7.5

BI8

0

0.02

0.05

2.5

BI9

0

0.01

0.02

2

BI10

0

0.05

0.1

2

Bioprinted tissues consolidation
Bioprinted tissues were consolidated for 90 minutes at 37°C in a solution containing at least one of the following components: CaCl2 (Sigma, France), transglutaminase (Ajinomoto Activa WM, Japan) and thrombin (Sigma, France). Table 2presents the composition of the different consolidation solutions used during the study. Once the consolidation completed, each bioprinted tissue was rinsed three times with sterile NaCl 0.9% (Versol, France) for 5 minutes at room temperature.
Table 2
Final composition of the 24 consolidation solutions used in the present study.
CODE

Ca2+ (mM)

Transglutaminase (mg/mL)

Thrombin (U/mL)

C1

270

40

0

C2

270

40

10

C3

27

40

0

C4

270

4

10

C5

270

80

0

C6

270

4.10− 2

0

C7

270

160

0

C8

270

4.10− 1

10

C9

1360

40

0

C10

2.7

40

0

C11

2700

40

0

C12

0

40

0

C13

0.27

40

0

C14

27

40

10

C15

0

40

0

C16

270

0

0

C17

1.8

40

10

C18

270

2

10

C19

27

4.10− 1

10

C20

27

4

10

C21

270

5.10− 1

10

C22

0

4.10− 1

10

C23

270

20

10

C24

2.7

20

10

Bioprinted tissue culture and monitoring
After the consolidation step, the different bioprinted tissues were grown in nontreated 6 or 12well plates at 37°C in a 5% CO2 atmosphere, in their dedicated culture medium (Supplementary Table S1). Spent media were recovered in triplicate three times a week and placed at 20°C until quantification of the lactic acid secreted by cell metabolism.
The lactic acid concentration was quantified using the LLactic Acid Assay from Megazyme (LLactic Assay Kit, KLATE, Megazyme, France) for autoanalyzer procedures. The assay was performed according to the manufacturer's instructions. Each sample was deposited in triplicate in the 96well plate. For highly concentrated samples, dilutions were performed in deionized water and respective commercial fresh medium was used as a negative control. The plates were incubated at 25°C in the dark for 10 min before measurement of the absorbance at 340 nm with a spectrophotometer (TECAN infinite®).
DMA measurements
The viscoelastic behaviour of the different consolidated hydrogels was characterised by frequency sweep experiments in DMA in compression mode at 37°C. These experiments were conducted with a rotational rheometer (DHR2, TA Instruments, Guyancourt, France) with a DMA mode (torque = 0N) using diskshaped samples and a parallel plate geometry (8 mm). A preliminary study was performed to define the linear viscoelastic domain, which corresponds to the displacement range where the material properties are assumed to be constant. This domain is determined using oscillatory compression experiments with constant frequency and varying displacement. Then, dynamic compression tests were performed with a frequency range of 0.1 to 10 Hz (i.e. 0.628 to 62.8 rad/s) at a constant displacement, which is within the linear viscoelastic regime. In these dynamic compression tests, bioprinted tissue undergoes a periodical mechanical strain \(\epsilon\) of less than 0.2%, \({\epsilon }_{0}\) and of angular frequency \(\omega\) following the Eq. (1):
\(\epsilon = {\epsilon }_{0}sinsin \left(\omega t\right)\)

Equation (1)

In the case of the generalised Maxwell model, the storage E′(ω) part of the complex modulus is expressed by the equations (2):
\({E}^{{\prime }}\left(\omega \right)= {E}_{0}\left[1+{\sum }_{\alpha =1}^{m}\frac{{\beta }_{\alpha }\omega ²{\tau }_{\alpha }^{²}}{1+\omega ²{\tau }_{\alpha }^{²}}\right]\)

Equation (2)

\(\eta = {E}_{0}{\sum }_{\alpha =1}^{m}{\beta }_{\alpha }{\tau }_{\alpha }\)

Equation (3)

where \({E}_{0}\) is the Young’s modulus of the isolated spring. The relaxation times, \({\tau }_{\alpha }^{}\), and the dimensionless reference parameters \({\beta }_{\alpha }\) stand for the contribution of each branch to the global modulus. The overall viscosity \(\eta\) can be defined as Eq. (3). The timeconstant values were regularly distributed between the reciprocals of the highest (62.8 rad/s) and the lowest (0.628 rad/s) angular frequencies of the experimental dynamic modulus. The chosen number of modes was sufficiently high to obtain accurate fitting, but not too large to avoid inconsistent results (e.g., negative values of \({\beta }_{\alpha }\)). Practically speaking, this led to threetime constants (m = 2), regularly spaced on a logarithmic scale between 5 × 10− 2 s and 5 × 10− 1 s.
Identification was achieved by solving the following minimization problem described by Eq. (4):
$${f}_{obj}^{}\left({E}_{0},{\beta }_{1}\dots {\beta }_{m}\right)={\sum }_{i}^{k}\left[\frac{({E{\prime }}_{i}{E{\prime }}_{i}^{mes})²}{{E{\prime }}_{i}^{mes}}\right]$$
4
where \({E{\prime }}_{i}^{mes}\) is the storage modulus obtained from the measured data and \({E{\prime }}_{i}\) is the one computed with viscoelastic parameters. \(k\) is the number of measurements acquired during the frequency sweep compression test. The optimization procedure was performed by using the Microsoft Excel Solver (version 2016) with the Generalised Reduced Gradient (GRG) nonlinear solving method.