Allometric Scaling of physiologically-relevant organoids

The functional and structural resemblance of organoids to mammalian organs suggests that they might follow the same allometric scaling rules. However, despite their remarkable likeness to downscaled organs, non-luminal organoids are often reported to possess necrotic cores due to oxygen diffusion limits. To assess their potential as physiologically relevant in vitro models, we determined the range of organoid masses in which quarter power scaling as well as a minimum threshold oxygen concentration is maintained. Using data on brain organoids as a reference, computational models were developed to estimate oxygen consumption and diffusion at different stages of growth. The results show that mature brain (or other non-luminal) organoids generated using current protocols must lie within a narrow range of masses to maintain both quarter power scaling and viable cores. However, micro-fluidic oxygen delivery methods could be designed to widen this range, ensuring a minimum viable oxygen threshold throughout the constructs and mass dependent metabolic scaling. The results provide new insights into the significance of the allometric exponent in systems without a resource-supplying network and may be used to guide the design of more predictive and physiologically relevant in vitro models, providing an effective alternative to animals in research.

Recent advances in in vitro technology have led to the development of organoids, i.e. cell aggregates grown from a small number of stem cells able to self-assemble, recapitulating the three-dimensional architecture of an organ at the micro-scale 1,2 . They are currently considered as one of the most promising ways to study cell behaviour and may have significant potential as a tool for investigating developmental biology, disease pathology, regenerative medicine and drug toxicity. Moreover, their downscaled functional and structural resemblance to mammalian organs suggests that they might also follow allometric scaling rules, such as Kleiber's law (KL) [3][4][5] .
KL is one of the most well-known allometric relationships in biology. It states that the basal metabolic rate (BMR) of an organism scales with its body mass, M, according to a universal power-law BMR = aM α , (where α is a scaling exponent, a a constant). On the basis that the number of cells in an organism is proportional to its mass, an analogous expression relates the average resource consumption rate per cell (cellular metabolic rate, CMR) to mass as CMR∝ M (α−1) = M b . The majority of experimental studies support the claim that KL applies across species with α = 3/4, or b = −1/4 hence it is often referred to as the quarter power scaling law 6 . A number of theoretical explanations have been proposed to support quarter power scaling 5,7 . As KL applies across several orders of magnitude of mass, it has been described as a unifying equation of biology 8 . Indeed, metabolic requirements of organisms, as well as their constituent cells, influence many fundamental biological properties common to all forms of life such as growth, differentiation and death. The pertinence of KL to the development of physiologically relevant in vitro models has been highlighted 9,10 , and it is considered by some scientists as the key to enabling the extrapolation of biological parameters from micro-scaled in vitro systems to the in vivo context: one of the holy grails of modern biomedical science 11 .
Ahluwalia 12 recently demonstrated that cells in vitro can maintain KL even the in absence of a resource-supplying network. In particular, given that oxygen is the primary metabolic resource for cells, the author evaluated the CMR of cell spheroids with increasing radii using computational mass transfer models coupling Michaelis-Menten (MM) reaction kinetics and oxygen diffusion through 3D constructs. The results show that KL holds in 3D spheroids if they have a suitable mass and cell density to enable the formation of a sufficiently large concentration gradient within the construct. However, the model proposed by Ahluwalia did not focus on organoid technology as it did not take into account the fact that organoids generated from embryoid bodies rely on the use of a gelatinous protein mixture resembling the extracellular environment surrounding the construct (e.g., Matrigel or Geltrex). While the 3D extra-cellular matrix (ECM) mimic supports cell survival and function [13][14][15] , it might be a limiting factor for oxygen uptake within the structure, particularly in the case of non-luminal in vitro replicates such as brain, liver or pancreatic organoids. Moreover, the cell density of standard in vitro cultures considered in the paper (i.e. 5.14•10 12 cells/m 3 ) is about two orders of magnitude lower than those typically reported for such organoids (around 10 14 cells/m 3 ) which are close to physiological cell densities 1,13,14,16 . Thus, for a fixed mass or radius, the oxygen consumption rate per unit volume of non-luminal organoids is much higher than in spheroids, which often leads to a dead core -a recurring feature in these systems 17 . By combining computational models with images of mature brain organoid slices stained with a nuclear dye, Berger et al. 17 recently identified a minimum threshold oxygen concentration for viability, C crit , of 0.04 mol/m 3 , which is also the critical value needed for mitochondrial ATP production 18 .
The scientific community is aware of the urgent need to improve organoid technology in terms of quality and reproducibility of sample size and shape. Particular attention should be addressed to eliminate, or at least minimise, the 'non-viable' centres while maintaining the ability of stem cells to self-organize and differentiate recapitulating physiological forms and functions [19][20][21][22] . This study is thus aimed at identifying a working window in which two criteria necessary for physiological relevance in cell-filled (non-luminal) organoids are met: (i) they follow "physiological" quarter power scaling; and (ii) the oxygen concentration within their volume is maintained above the critical threshold C crit . By customising the computational framework based on allometric scaling outlined by Ahluwalia 12 , the CMR of spherical organoid models was computed over 5 orders of magnitude of mass to estimate the value of the exponent b. Additionally, the concentration gradient through the models was mapped to determine their minimum oxygen levels. Finally, to investigate the effects of integrating vasculature in organoids, a single fluidic channel was added to the models.

computational Methods
To compute the oxygen concentration in organoid systems and determine the allometric relationships between the cellular metabolic rate (CMR) and organoid mass, 3D multi-physics models were implemented and solved numerically using COMSOL Multiphysics (version 4.3 COMSOL AB, Stockholm, Sweden) and the UMFPACK direct solver. In the models, reaction, convection and diffusion of oxygen through organoids and gels as well as fluid dynamics for incompressible fluids were considered. organoid volume expansion and stem cell proliferation. Referring to published data on brain organoids derived from human induced pluripotent stem cells (hiPSCs) 13,17 , we modelled both stem cell proliferation and organoid volume expansion to identify a range of experimental organoid masses and corresponding cell densities. In particular, cell proliferation was modelled using the Sherley model 23 : where P(t) is the population size at time t, P 0 the initial number of cells, γ the mitotic fraction and DT is the division time. Since we supposed that all the cells within the organoid reach mitosis, γ was set to 1. The division time, DT = 1.19•10 5 s, was estimated from the data reported in Berger et al. 17 , where 9000 hiPSCs were plated and a total of 1.8•10 5 cells were counted at the end of the experiment after 20 days. On the other hand, the function describing organoid expansion was extracted by least squares fitting the data on measurements of organoid diameter reported Monzel et al. 13 to a generic exponential: where r(t) is the organoid radius at time t, r min the initial radius (i.e. 250 μm) and r max the maximum radius reached at the end of its growth (i.e. 562 μm). The expansion coefficient, λ, was estimated as 1.126•10 −5 s −1 (R 2 : 0.9884, RMSE: 16.43). From the equations above, we calculated the cell densities at three time points (i.e. day 4, day 8 and day 12, Table 1). Since the cell densities vary by less than 20% over time, a mean value of 2.52•10 14 cells/m 3 was used for all three time points. This value is of the same order of magnitude as those reported in other protocols for non-luminal organoid generation (e.g., ref. 24 ).
When developing the computational models, four basic phenomena were considered: (i) cell proliferation, (ii) organoid volume expansion, (iii) oxygen consumption within the organoid and (iv) oxygen diffusion. Given that the characteristic times for organoid growth and expansion are much longer than the times for oxygen diffusion Organoid density = 1000 kg/m 3 . Data on hiPSC generated brain organoids from Monzel et al. 13 and Berger et al. 17 .
www.nature.com/scientificreports www.nature.com/scientificreports/ and consumption within the spheroid (Table 2) Table 1  www.nature.com/scientificreports www.nature.com/scientificreports/ Cerebral organoids are not perfect spheres, but slightly ellipsoidal. In Fig. 1A,B, the spherical constructs can be considered as the worst case in terms of surface area to volume ratio -and hence minimum oxygen concentration. On the other hand, the model in Fig. 1C represents a structure with a blood vessel mimic, similar to what might occur in implanted cortical organoids after they connect to a host vasculature 21 . the multi-physics models. A reaction and diffusion model 25 Table 1), we assumed that the oxygen diffusion constant within the organoid (D org ) is the same as that in in vivo tissues 27 . Oxygen consumption was described using MM kinetics, considering literature values of k m and maximal oxygen consumption rate (OCR) typical for human stem cells (see Table 3).
The boundary oxygen concentration was fixed at 0.2 mol/m 3 , assuming a well-mixed supply of medium to the organoids. The ruling equation for this system is: More detailed equations and boundary conditions are reported in the Supplementary Materials. In the second configuration in Fig. 1B, only diffusion was implemented within the hydrogel shell. The oxygen diffusion constant in the shell, D M , was equal to that measured for Matrigel 28 . The cell-filled core was identical to the basic model in Fig. 1A.
Finally, in the third configuration, oxygen mass transport considering diffusion, convection and consumption in the spheres (see Supplementary Materials) were coupled to fluid dynamics in the central channel based on the geometry shown in Fig. 1C. The velocity field in the channel was solved using the Navier-Stokes equation for Newtonian fluids, setting the inlet velocity to 10 −3 m/s, a value consistent with blood velocity within capillaries 29 . The simulations were run on constructs with the same radii as in ref. 12 (ranging from 54.5 to 2142 µm) and includes the experimental radii reported in Table 1. All the constants and variable ranges used are listed in Table 3.
Data evaluation and fitting. Once the solutions were obtained, the minimum oxygen concentration within the organoid domain was computed. In addition, the surface integration function on the organoid domain was used to determine the total inward oxygen flux at the organoid boundaries and the resulting CMR as in Ahluwalia 12 . Specifically, the average metabolic rate per cell or CMR is given by the surface integral of the inward oxygen flux, divided by the number of cells in the construct. Both the CMR data and the minimum oxygen concentrations were then imported into Matlab 2017a (The Mathworks Inc., Boston Massachusetts) for plotting and curve fitting. To investigate the allometric behaviour of the constructs, log(CMR) versus log(mass) graphs were generated for each configuration using the method described in ref. 12 and outlined in the Supplementary Materials. The allometric exponent b was estimated from the slope of the log-log curves for the range of experimental masses in Table 1. The best fit value of b, corresponding to the slope covering the maximum number of points and with the maximum R 2 value, was also determined. Here, the log(CMR) -log(mass) dataset for each configuration was fitted to a straight line discarding one point at  Table 3. Parameter inputs to the computational models.
a time beginning from the lowest mass. A Gaussian equation was used to fit the minimum oxygen concentration data for each configuration implemented. The equations were used to estimate the radius at which the oxygen concentration is equal to C crit in the organoid domain.

Results
As shown in Fig. 1A, for gel-free models with small masses, the CMR is mass-independent (i.e. the slope b = 0). This implies that all the cells in the constructs consume metabolic resources at the same rate since they all perceive the same oxygen concentration 12 . In the range of experimental organoid masses reported in Table 1, the slope b = −0.06 (R 2 : 0.95), which is far from quarter power scaling. Larger masses (corresponding to radii between 600 and 2000 µm) have a slope close to the KL value, b = −0.23 (R 2 : 0.99, best fit). The model suggests that a resource-supplying network is not a necessary condition for quarter power scaling, but allometric scaling will not occur below a minimum construct mass. The constructs surrounded by a shell also have a constant CMR for small masses, although its value is slightly lower than that of the gel-free ones (Fig. 1B versus 1A). This is reasonable because the oxygen concentration at the surface of the organoid domain with a gel shell is always less than 0.2 mol/m 3,12 . In the experimental range of interest b is close to −1/4 as in KL (b = −0.23, R 2 : 0.98), falling to b = −0.33 (R 2 : 0.99, best fit) for larger spheres. In the case of the constructs with the hydrogel shell and the micro-channel, the log-log graph for CMR versus mass (Fig. 1C) shows that the experimental range of masses also follow KL as the slope b is almost −1/4. Moreover, b is still mantained within the allometric range for larger masses (b = −0.25, R 2 : 0.98). Figure 2 reports the minimum oxygen concentration values, C min , within the organoid domains as a function of their radius and the fitted Gaussian curves for the three configurations represented in Fig. 1. The viability threshold oxygen level in brain organoids C crit 17 is also indicated in the graph. For a given radius, the minimum oxygen concentration in the gel-free constructs (blue) is higher than for the hydrogel-encapsulated ones (red). The radii at which the minimum oxygen concentration is equal to C crit determines the size limit for these systems, and corresponds to 712 µm for the gel-free models and 475 µm for the organoids surrounded by gel. On the other hand, comparing the green and red curves, it is clear that the use of a microfluidic channel can improve oxygenation within the constructs even in the presence of an ECM-mimicking shell. In fact, the maximum construct radius in the presence of a microchannel is 616 µm. Table 4 summarizes the range of radii for which both criteria, b≈−1/4 and C min ≥ C crit are satisfied for each configuration presented in Fig. 1. From the table it is clear that the working window for physiological relevance is very narrow for the organoid models sheathed by an ECM gel. Here, only constructs with radii between 236 and 475 µm satisfy both criteria. For the gel-free models, the range is even narrower, going from 654 to 712 µm. A much larger range of radii can be covered if a single microfluidic channel is incorporated in the organoids, as both b≈−1/4 and C min ≥ C crit are valid for the range 236 to 616 µm.

Discussion
Developing physiologically relevant in vitro models to mimic cell and tissue organ response is one of the main challenges in experimental cell biology and tissue engineering. Organoid technology has emerged as a novel method to generate structural and functional tissue models 15,[30][31][32][33] . Luminal organoids were first generated from intestinal stem cells by Clevers and coworkers 34 . As they are hollow, they are unlikely to suffer from necrotic cores 35 . On the other hand, non-luminal organoids such as brain and liver often possess non-viable centres, principally due to the lack of oxygen which is widely considered as the limiting metabolic resource for living cells 17,22 . Local oxygen levels also influence cell fate [36][37][38][39][40] , likely compromising tissue function and structural arrangement. Thus, further strategies to guarantee the nutritional needs of such organoids are required.
It has been suggested that allometric behaviour could be considered as a benchmark of physiological relevance 9,25 and a recent study shows that KL can be achieved in high density constructs 12 . The characteristic features www.nature.com/scientificreports www.nature.com/scientificreports/ of non-luminal organoids, that is high cell densities and their potential to therefore follow physiological metabolic scaling as well as to adopt the structural and functional features of mammalian organs on one hand, and their tendency to develop an oxygen deprived core on the other, appear to be in contrast. In this context, our aim was to evaluate the minimum oxygen concentration and the CMR of non-luminal organoids to enable the design of cell culture protocols in which both quarter power scaling holds and the oxygen level in the organoids is greater or equal than the necrotic threshold of 0.04 mol/m 3 (C crit ). As illustrated in Fig. 3, the intersection of the two criteria is argued to represent the window of physiologically relevant organoid sizes for designing constructs with better translational value. We therefore established 3D models of cell-filled spheres in a variety of configurations covering over 5 orders of magnitude of mass and including the range of sizes reported in the literature. The cell density and oxygen consumption parameters were those reported in the literature for brain organoids and stem cells respectively, while the boundary oxygen concentration was set at 0.2 mol/m 3 , as in ref. 12 . To identify the range of masses in which KL holds, the slope of log-log graphs of CMR versus mass was estimated. The masses and corresponding radii for which all the cells in the 3D constructs receive an adequate supply of oxygen were also computed.
Some mini-organs, typically liver buds generated using mesenchymal cell-driven condensation, are dependent on an ECM-like substrate but do not require embedding in a gel to undergo differentiation and self-organisation 24 .
Our results indicate that if cultured in the presence of a well-mixed and constant supply of freshly dissolved oxygen at the organoid-media boundary, these structures do not suffer from nutrient limitations over a wide range of masses (up to 1.5 mg, Table 4). The results are confirmed by experimental studies on liver buds by Ramachandran et al. 22 ; the authors found no evidence of cell death in liver buds cultured in well-mixed fluidic bioreactors while those in static conditions developed a necrotic core. Mattei et al. 33 demonstrated that the minimum oxygen level in liver buds with an average mass of a few mg cultured in bioreactors is 0.04 mol/m 3 , compared with 0.01 mol/m 3 in static culture plates. As shown in Table 4, the gel-free structures resembling liver buds maintain the criteria of quarter power scaling (b ≈ −1/4) and viable oxygen concentrations only if their radii are maintained in a narrow range (between 654 and 712 µm in the models).
On the other hand, most of the current protocols for generating brain organoids rely on the use of ECM gels (e.g. Matrigel) which act as a 3D media supporting cell proliferation and differentiation. However, the gels also limit oxygen diffusion within the sample: indeed non-viable cores have been reported in these organoids. In fact, in the models the radii at which both the conditions (i.e. allometric scaling and minimum oxygen concentration above C crit ) are satisfied are smaller than those for the gel-free constructs. Interestingly, in most studies, the sizes of organoids surrounded by a shell of ECM are generally smaller than those without a shell 13,14,22,33 .
We also developed a model to determine if the incorporation of a simple fluidic channel mimicking an elemental central blood vessel could widen the working window in which both criteria of b = −1/4 and C min ≥ C crit are satisfied (Fig. 1C). The results show that the addition of the channel results in oxygen concentrations above C crit , and quarter power scaling for a wider range of radii, including those for brain organoids summarised in Table 1. Such channels have been reported in cortical organoids which are integrated into the host vasculature when implanted in animals 21 , but their fabrication in vitro remains a challenge.  Table 4. Range of construct radii in which both KL law and oxygen concentration above the 0.04 mol/m 3 threshold are maintained.

Figure 3.
Physiologically relevant organoids should obey allometric scaling laws and cells in the centre of the constructs must be exposed to oxygen levels above the critical threshold for vitality. Using computational models to identify the range of radii which respect these two conditions, we can establish criteria for generating in vitro models with better translational value.
It is important to highlight that we assumed a constant oxygen concentration at the outer surface of the geometries implemented, corresponding to the maximum dissolved oxygen concentration in water at 37 °C at atmospheric pressure: this condition requires well-mixed and continuously renewed media through the use of fluidic devices such as bioreactors. Furthermore, although current technology using oxygen-sensitive fluorescent probes have made it possible to determine local intra-cellular oxygen levels 35,41 , accurate measurements of single cell and cell aggregate resource consumption rates are difficult to perform and scarce in the literature. The maximum oxygen consumption rate and the MM constant used here are based on published data for human mesenchymal stem cells. They may vary slightly from one cell population to another as different cells have different metabolic demands, and the critical range of radii for physiologically relevant organoid systems may vary accordingly.
Based on the considerations in this work and our previous study 12 , typical values of b and their significance in terms of resource availability are reported in Table 5. The exponent b is close to zero when all cells in the volume consume resources at the same rate because they all perceive the same oxygen concentration. As the constructs increase in size, cells on the periphery perceive and consume higher concentrations of oxygen than cells close to the centre, which adapt to lower levels of oxygen and hence lower consumptions rates -a characteristic of MM kinetics. As a consequence, the average CMR decreases and at the 'knee' of the log-log CMR versus mass graph, b is approximately −1/4 as in KL. At higher masses, as shown in ref. 12 , the coefficient b tends to −1/3, which indicates that resource uptake depends on the surface area of the construct 42 and is diffusion limited.
In conclusion, we suggest that a resource-supplying network is not a necessary condition for quarter power scaling in high density organoids or cellular aggregates. However, despite obeying KL, non-luminal organoids may be resource limited as core oxygen levels fall below C crit . Thus, quarter power scaling in these structures is not a sufficient criterion for physiological relevance because part of the cell population may not perceive enough oxygen to guarantee their vitality. The range of sizes which guarantee physiological relevance is in fact quite limited and should be taken into consideration when employing organoids as mini in vitro models of human organs.
There is a great demand for improving organoid models with respect to their size, structure and composition: this study provides a first step towards the development of protocols for designing more predictive organoid systems based on objective criteria. Although current methods for integrating organoid and fluidic technology have focused on on-chip devices in which constructs are cultured inside fluidic channels 24 , we suggest that efforts be made to develop methods using bioprinting technology to fabricate organoids with a perfusable central 43 . The design of simple microcirculation-on-a-chip solutions must however be accompanied by accurate methods for regulating and monitoring oxygen and its consumption rates.

Data Availability
The equations and tables provided in the manuscript are sufficient to replicate the computations. Data files are available upon reasonable request.   Table 5. Typical values of the allometric exponents and their significance in terms of resource availability. (N.B. α = b + 1).