Quantitative characterization of 3D bioprinted structural elements under cell generated forces

With improving biofabrication technology, 3D bioprinted constructs increasingly resemble real tissues. However, the fundamental principles describing how cell-generated forces within these constructs drive deformations, mechanical instabilities, and structural failures have not been established, even for basic biofabricated building blocks. Here we investigate mechanical behaviours of 3D printed microbeams made from living cells and extracellular matrix, bioprinting these simple structural elements into a 3D culture medium made from packed microgels, creating a mechanically controlled environment that allows the beams to evolve under cell-generated forces. By varying the properties of the beams and the surrounding microgel medium, we explore the mechanical behaviours exhibited by these structures. We observe buckling, axial contraction, failure, and total static stability, and we develop mechanical models of cell-ECM microbeam mechanics. We envision these models and their generalizations to other fundamental 3D shapes to facilitate the predictable design of biofabricated structures using simple building blocks in the future.

show the differences between 3D printed and cast collagen networks below in Supplementary

Buckling Wavelength Measurements
Qualitative observations of the cell-ECM microbeams reveal that beams originally printed straight develop undulations over a 24-hour period. To quantify the wavelength of the undulations and test EB beam theory, we take 3D confocal fluorescence microscopy images of the cell-ECM microbeams, using fluorescently labelled cells (see Methods). We employ three different methods of analysis to measure the wavelengths from the 3D confocal stacks: intensity autocorrelation analysis, 3D curve fitting of beam backbones, and direct visual inspection. We find consistent results with all three methods.
To determine λ when the undulation amplitude is small, we rely on an intensity autocorrelation analysis that is sensitive to regular oscillations. First, each confocal stack is filtered to remove noise using a 3D Gaussian kernel, one pixel in width. The filtered stack is rotated in 3D using image moments analysis to ensure that the long axis of the beam is parallel to the closest Cartesian axis, which we define as z. The filtered and rotated stack is then digitally stretched or compressed to ensure cube-shaped voxels, creating an intensity distribution, I(x,y,z), where x and y correspond to the axes transverse to the beam direction. A maximum intensity projection is then taken along the y-axis which corresponds roughly to the direction of confocal scanning. From this projected image, we compute an intensity-intensity autocorrelation function using the Fourier and identify the location of the first peak at finite Z, which equals the undulation wavelength, λ ( Supplementary Fig. 3c).
Supplementary Figure 3. (a) A 2D projection of a confocal fluorescence 3D stack shows the shape of a buckled microbeam with a small undulation amplitude. We show a digitally stretched image to accentuate the undulations (spline drawn manually to guide the eye). (b) To measure the buckling wavelength, λ, we compute a 2D intensityintensity autocorrelation function, C(X, Z). The 2D correlation function shows peaks located at Z = ±λ. (Scale bar: 100 µm) (c) A plot of C(0, Z) is used to measure λ, where the location of the first peak at finite Z is equal to λ.
To measure the undulation wavelengths using a different method for comparison, we preprocess the confocal stacks in the same way as described above, but instead of computing a correlation function, the backbone locations of the beam are measured in the x-y plane while stepping through the beam's long axis, z. The x-y centroids at each z location are found by computing the 2D centre of intensity 4 . Prior to computing the centroids, we fill in the void-space between the cells by carefully dilating and eroding the images producing a closed continuous cross-section. This step of processing is performed by direct visual inspection for each beam to carefully choose appropriate thresholding, dilation, and erosion parameters for each beam. The resulting backbone locations are fit simultaneously to two sinusoidal functions, given by We initiate the fitting procedure by visually matching the observable undulation wavelength and amplitude, λ1 and A1 to the backbone data; the second set of parameters, λ2 and A2 account for any weak bends corresponding to a λ2 much larger than the length of the beam. While this procedure requires a lot of manual intervention, it produces measurements of λ1 that agree with the correlation analysis described above as well as the undulations that can be observed by eye ( Supplementary   Fig. 4). Whether using correlation analysis or curve-fitting, all wavelengths are compared to those measured by direct visual inspection by identifying the locations of peaks and troughs in transverse-projections of the confocal fluorescence stacks; all results agree to within the spread of the data shown in Figure 3d in the manuscript body.
Supplementary Figure 4. The backbone locations of a buckled microbeam, measured by image centroiding (symbols), are fit by a sinusoidal function to determine the buckling wavelength. This method complements the correlation function analysis described above and produces results consistent with direct visual inspection.

Collagen Microbeam Elastic Moduli
To demonstrate control over collagen gelation, we perform low-amplitude oscillatory rheology on collagen gels cast between roughened 25 mm plates using an Anton-Paar MCR-702 rheometer. To mimic the solution conditions of our printing experiments, we prepare collagen-1 solutions in cell growth media supplemented with HEPES buffer (25 mM), and bring the pH to 7.4. The collagen solution is loaded into the rheometer geometry at 20 o C, and a small enclosure with a solvent trap is put in place to prevent evaporation during testing. A CO2 line is inserted into the enclosure to maintain 5% CO2, and the temperature is raised from 20 o C to 37 o C at a rate of 1 o C/min. After 30 minutes, we apply a 1% oscillatory strain at 1 Hz and measure the elastic response, G'. We find that the scaling between G' and collagen concentration, c, is consistent with a c 2.2 power-law for concentration between of 0.5 mg/mL and 2.5 mg/mL, which agrees with prior work done using similar curing protocols 5 . At the lowest collagen concentrations tested, this scaling relationship breaks down, yet we find that G' > G'', indicating the formation of system-spanning networks even at low concentrations (Fig. 3c).
To test whether the collagen network shear modulus, G', is the modulus relevant to the types of deformations observed in the 3D printed microbeams, we employ an alternative method to determine the elastic modulus; we 3D print cell-free collagen beams into the microgel media and manually apply axial loads while observing the response with photographic or microscopic imaging ( Fig. 3a, b). The macro-beams are deformed by pressing on one end with a microindentation system (Hysitron Bio-soft). We find that the macro-beams buckle and then straightenout again after the load is removed. Using the image analysis methods described above, we measure the bucking wavelength and beam-diameter; with these measurements combined with the known G' of the microgel medium, we are able to infer the beam elastic modulus, E. We extend this method to the micro-scale, performing a similar procedure on the microscope while using a printing needle to deform the microbeams (Fig. 3a, Supplementary Movie 2). Together, these tests on cell-free beams show that at the higher end of collagen concentrations, greater than 1.0 mg/mL, the elastic modulus relevant to buckling is very close to the collagen shear modulus. By contrast, at collagen concentrations of 1.0 mg/mL and below, we see that the datasets diverge ( Fig. 3c).
Since we are unable to measure E at concentrations below 0.5 mg/mL using this method, any experiments using collagen below this concentration are excluded from analysis if knowledge of E is required.
To investigate the differences between G' and E at low collagen concentrations, we image collagen networks at various concentrations, both cast on glass slides without surrounding microgel media and printed into the microgel medium. From laser-scanning confocal reflectance microscopy images, we find that collagen networks embedded in microgel medium look qualitatively similar to those cast as free liquid drops but have a visibly larger mesh-size at the lower concentrations of 0.25 mg/mL and 0.5 mg/mL. Since free collagen molecules can diffuse into the surrounding microgel medium during polymerization, it is possible that the competing rates of polymerization and diffusion favour diffusion at very low collagen concentrations, resulting in a sparser network than expected. While detailed investigations are needed to determine the underlying mechanism responsible for the differences between G' and E, our mechanical measurements and reflectance images produce consistent results. Moreover, all the analyses throughout this manuscript employ the measurements of E, avoiding the need for a theoretical model or detailed understanding of the polymerization process. The consistency in rheology and apparent network structure also indicates that the potential roles of viscoelasticity and network complexity may contribute comparably to both sets of measurements. We would expect viscous stresses to dominate at high shear rates; given that cell-driven beam buckling is a very slow process, we believe viscous stresses are likely small compared to elastic stresses or shear-rate independent yielding of microgels. Over long time scales the non-elastic contributions would likely arise from network remodeling and reorganization, which we have explored in multiple different control experiments described throughout the manuscript. We use only measured quantities to test the applicability of EB theory to cell-driven buckling; λ is predicted using measurements of E from cell-free collagen beams, combined with measurements of G', and R. Finding excellent agreement, we extend this approach to cell-ECM microbeams, using EB theory to determine E from measurements of λ. Using 57 measurements of λ from cell-ECM microbeams containing 3t3, GL261, and Panc02 cells, we compute the beam elastic modulus containing cells, E + , given by and plot the results versus the cooresponding beam elastic modulus without cells, E -, measured as described above (Supplementary Fig. 6). The 4 th power dependence of E on λ will lead to large variations in E + and Earising from small variations in λ, as observed in this plot. However, these datapoints span several decades in both directions, revealing a linear relationship given by E + = 1.13 E -. We use this relationship to determine E + in calculations of σcell where no λ is measured (break-up threshold and axial contraction datapoints), converting from Eto E + . In experiments where beams buckle and λ is measured, we directly solve for E + without using Ewhen calculating cell generated stress, σcell (see manuscript body describing models involving σcell).
Supplementary Figure 6. We 3D print microbeams having different collagen concentrations and measure their elastic moduli with and without embedded cells (E + and E -). In both cases, the buckling wavelength, λ, is used to determine the beam elastic modulus at a given collagen concentration. While the 4 th power dependence of E on λ generates large variations in E from small variations in λ, the large range of moduli measured shows a linear scaling, given by E + = 1.13 E -.
Recent investigations of collagen gel elasticity have shown these networks having strong non-linear responses to compression and stretch 6,7 . For example, the shear modulus increases strongly in stretched networks and decreases weakly in compressed networks. Correspondingly, the modulus of uniaxial strain increases strongly with stretch along the same axis and decreases weakly with compression. One major outcome of these observations is the discovery of a highly asymmetric Poisson's ratio about the unstrained state. While these strong nonlinearities have major implications in the mechanical performance of collagen networks in vitro and in vivo, it is unclear how they may contribute to the mechanics of the cell-ECM microbeams investigated here. For example, when cells within the microbeams contract they probably initially tense the collagen networks, suggesting that stiffening may occur. However, in all cases in which microbeams buckle, fail, or contract, the net effect is to contract the network; the pliable microgel environment allows the entire collagen network to compress, suggesting a weak strain softening could be occurring.
We hope that our results inspire new investigations of collagen elasticity in contexts where the network boundary conditions evolve in time. We believe that using jammed microgels as encasing environments for biopolymer networks represents a path forward for such studies with tuneable boundary conditions.

Wavelength measurements broken down by cell type and test type
We compare various measurements of λ to predictions from EB theory, using the classic result, λ = 2π (EI/G') 1/4 . Here we determine E, I, and G', from experimental measurements described above and in the manuscript body. For the three cell types investigated here, and for cell-free beams under manual loading, we find that EB theory predicts our data with no fitting parameters (R 2 = 0.93, Supplementary Fig. 7). Control experiments without cells exhibit no spontaneous buckling, indicating that undulations emerging in cell-ECM microbeams are a form of buckling driven by the contractile stresses generated by the cells within.

Elastic Shear Modulus and Yield Stress of Microgel Medium
To characterize the material properties and yielding behaviour of the mircrogel media used here, we perform both oscillatory and unidirectional shear rheology on samples prepared at several microgel concentrations (2.2% to 10%) using a Malvern Kinexus rheometer equipped with a roughened cone on plate geometry. Within the oscillatory frequency range of 10 -3 -10 2 Hz, we find that G' is greater than G'' for all polymer concentrations, indicating that the samples are dominantly solid-like ( Supplementary Fig. 8a, b). Choosing 1 Hz as a representative frequency, we plot G' versus microgel concentration, c, finding a scaling relationship consistent with G' ~ c 9/4 , as previously reported for the same microgels prepared in water 8 . To measure the yield stress of these microgel samples, we perform unidirectional shear tests described previously 1,2,8,9 . We find the experimental relationship between stress and shear-rate is well described by the Hershel- (1 ( / ) ) σ σ γ γ = +   , where σ is applied stress, σy is yield stress, γ is shear rate, c γ is the crossover shear-rate between solid-like and fluid-like behaviours, and p is a dimensionless number of order 0.5 ( Supplementary Fig. 8c) 10 . We determine σy for each sample by fitting the Hershel-Bulkley model to the data, finding the scaling relationship between σy and c consistent with σy ~ c 9/4 , just as found with the G' data ( Supplementary Fig. 8d). Plotting σy versus G' reveals an excellent correlation between the two variables, in very good agreement with our previous findings that σy = 0.13 G', which we found to arise from the elastic deformation energy required to drive single-particle rearrangements during shearing (Supplementary Fig. 8e).

Microscopy at the microgel-collagen interface
To examine the interface between collagen microbeams and the surrounding microgel medium, we print 1.9 mg/mL collagen-1 solutions into 0.15% (w/w) microgel medium. To image the microgel medium, we disperse 4 µm diameter fluorospheres which act as fiducial markers.
Imaging with laser scanning confocal fluorescence microscopy, we examine y-z cross-sections by integrating the intensity from a 3D stack along the beam's long axis (the x-axis). From these y-z projections, we see a circular void that shows where the microgels have been displaced by the deposited collagen. To quantify the sharpness of the interface between collagen and microgels, we examine the x-y slice through the beam's mid-plane. By integrating the intensity of this x-y slice along the x-axis, we see how particle density rises from within the interface and plateaus in the bulk, far from the beam. This intensity profile rises from 0 to its plateau value over a distance of about 25 µm. Since these microgels are approximately 4-5 µm in diameter, we expect intermixing between collagen and several layers of microgels to occur within this zone.

Length-dependence of microbeam contraction
To test whether friction may dominate beam contraction in a long-beam limit, we print cell-ECM beams of different lengths (1.0 mg/mL collagen) into microgel media with a yield stress of 1.9 Pa.
This combination of parameters lays in the contraction regime displayed in Figure 4b. Like the 10 mm long beams described in the manuscript body, we find that short beams, 1 mm in length, also contract. By contrast, much longer beams, 30 mm in length, do not contract at all. We find that the beam appears to compact radially and becomes smoother at its surface over a 24 hour delay between measurements, yet identifiable features do not appear to move at all along the axial direction. Thus, a friction-dominated limit appears to emerge with increasing beam length, which we hope to explore further in future work ( Supplementary Fig. 10). plots of C(R) × R 1.95 exhibit oscillations at longer length-scales that arise from cell-cell spatial correlations ( Supplementary Fig. 9b). By measuring the location of the first peak in these plots, we determine the average nearest-neighbour distance, Rnn, between cells in the microbeams. For the lower volume fraction beams, with φ < 5%, we can validate this procedure by visually counting the number of cells and measuring the dimensions of beams with ImageJ, finding agreement in volume fraction to within a factor of two. Additionally, Rnn measurements were checked visually using ImageJ by identifying nearest neighbours and measuring their centre-to-centre distances by eye, finding agreement to within 31%. These errors in measurements of linear dimension and volume fraction are consistent with one another; an error of nearly 30% in a linear dimension will lead to roughly a factor of two error in volume approximations. While a higher degree of certainty is desirable, since the data span multiple decades in parameter space, factors of two in either direction do not change the principle conclusions.

Supplementary
Supplementary Figure 11. (a) The nearest neighbour distance between cell pairs, Rnn, is measured by computing the intensity-intensity autocorrelation function, C(R), of confocal stack intensity distributions. We find that C(R) resembles a modified Lorentzian distribution decaying to a power of 1.95. (b) Multiplying C(R) by R 1.95 removes the decay and accentuates long-wavelength oscillations that arise from cell-cell spatial correlations. We use the location of the first peak in this plot to determine Rnn.

Literature data points
The datapoints shown in Figure 5b labelled "2D literature" were taken from published manuscripts describing traction force microscopy measurements of cells on soft, flat, 2D substrates [11][12][13]  reported, respectively. We treat these measurements as the 2D equivalent of the cell generated stress estimated from our data and plot them without modification onto our graph of cell-generated stress in 3D (Fig. 6b). We find that a scaling law fit to our data in 3D extrapolates through this 2D traction stress data to within a factor of approximately 2. Additionally, recent investigations of cell-generated traction forces in 3D matrices agree very well with our measurements, laying close to the extrapolated fit to our data-points [14][15][16] . From the traction forces distributions and strain energies measured using 3D TFM methods in the literature, we estimate the cell-generated stress ( Supplementary Fig. 12 Figure 12. We surveyed the literature for measurements of traction-stresses exerted by cells on 2D elastic substrates. A scaling law fit to our 3D data extrapolates to the 2D data when rescaled by a factor of 2.25. Threedimensional traction force microscopy literature were also surveyed. Shown in open symbols, we establish that these data also fall on our prediction of cell generated stress relative to microenvironmental modulus.

3D Printed Neural Crest Model
To test our 3D bioprinting method's ability to generate a model developing tissue made from multiple cell types or materials, we designed a 3D analogue of textbook-level illustrations of the developing neural crest. Fetal neural development is typically depicted as a series of snapshots showing the formation of the neural tube through the buckling of the neural plate and the translocation of the neural plate border and epidermis. In such a time-series, the epidermal sheet converges to cover the neural tube. To demonstrate our method's capabilities we 3D print structures made from fluorospheres (Fig. 7). The same structure can be printed from living cells, though we recognize that mimicking morphogenesis in vitro using such structures is beyond our current capabilities. To fabricate the 3D models, we fill three separate wells with solutions of fluorospheres of different colours. We program our printing stage to visit each separate well and draw an arbitrary volume of fluorosphere solution into a syringe. The different fluorosphere solutions are printed one at a time; when changing from one colour to the next, the syringe is programmed to empty into a waste well and wash itself with clean water and then clean microgels.
This method allows the printing of one seamless structure out of three different materials with a single set of printing commands. The green beads represent the epidermis, the blue beads represent the neural plate border, and the red beads represent the neural plate. We envision that the ability to fabricate these structures from living cells will one day enable testing a watchmaker's approach to tissue engineering, in which structures at each stage of development are created and allowed to further mature, with the goal of discovering how to set in motion later stages of development from conditions mimicking earlier stages.