## Abstract

Oxygen plays a central role in cellular metabolism, in both healthy and tumour tissue. The presence and concentration of molecular oxygen in tumours has a substantial effect on both radiotherapy response and tumour evolution, and as a result the oxygen micro-environment is an area of intense research interest. Multi-cellular tumour spheroids closely mimic real avascular tumours, and in particular they exhibit physiologically relevant heterogeneous oxygen distribution. This property has made them a vital part of *in vitro* experimentation. For ideal spheroids, their heterogeneous oxygen distributions can be predicted from theory, allowing determination of cellular oxygen consumption rate (OCR) and anoxic extent. However, experimental tumour spheroids often depart markedly from perfect sphericity. There has been little consideration of this reality. To date, the question of how far an ellipsoid can diverge from perfect sphericity before spherical assumptions break down remains unanswered. In this work, we derive equations governing oxygen distribution (and, more generally, nutrient and drug distribution) in both prolate and oblate tumour ellipsoids, and quantify the theoretical limits of the assumption that the spheroid is a perfect sphere. Results of this analysis yield new methods for quantifying OCR in ellipsoidal spheroids, and how this can be applied to markedly increase experimental throughput and quality.

## 1. Introduction

Oxygen plays a seminal role in cancer treatment and patient prognosis. The presence of molecular oxygen in a tumour markedly increases radio-sensitivity, with well-oxygenated regions responding to radiotherapy by up to a factor of 3 relative to anoxic sub-volumes [1,2]. This oxygen enhancement ratio is also seen in emerging modalities such as proton therapy [3,4], raising the tantalizing prospect of dose painting, where dose is selectively boosted to hypoxic regions to boost therapy response [5]. The basic idea underpinning dose painting has been discussed for over a decade, but application has been hampered by difficulty in non-invasive hypoxia imaging. Methods such as F-MISO PET (fluoromisonidazole positron emission tomography) have a maximum resolution in the millimetre regime, while oxygenation varies over a micrometre scale. As a consequence, mathematical modelling is vital for bridging the resolution gap [6].

Aside from therapeutic considerations, oxygen has a marked impact on patient prognosis. The pioneering work of Gray and colleagues in the early 1950s established that tumour oxygen concentration was correlated with prognosis, and extensive hypoxia was a negative prognostic marker [7]. This finding has been well replicated to present day [8–10] and is not solely due to hypoxia-induced treatment resistance. Under severe hypoxia, tumour cells can respond to such pressure by activating oxygen-sensitive signalling pathways [11,12]. Current biological thinking suggests that these signalling pathways act to alter gene expression to promote cell survival under adverse conditions. Hypoxia is also a major driver of angiogenesis, giving rise to new routes for cells to travel along [13,14], endowed with the ability to metastasize [15].

The extraordinary importance of oxygen in cancer treatment and evolution has made it an important avenue of study, with an urgent need for further research. Despite the fundamental importance of molecular oxygen in tumours, investigations have been complicated by the significant experimental difficulty in ascertaining oxygen concentration *in situ* [6]. Real tumours have highly heterogeneous oxygen supply and complex tortured vasculature, and even well-oxygenated regions are frequently inter-spaced with pockets of anoxia [14,16]. Standard two-dimensional monolayers of cells are not an ideal experimental model, typically exhibiting an unrealistically homogeneous oxygen contribution. There is however a more realistic experimental option in the form of tumour spheroids. These clusters of cancer cells grow in approximately spherical three-dimensional aggregates, and exhibit signalling and metabolic profiles more similar to real tumours than is observed in monolayer approaches [17–19].

Like monolayers, spheroids are relatively easy to culture, and growing interest has seen them used for a variety of purposes, including radiobiological application as a means to test fractionation [20–23], as a model for drug delivery [24–28], for investigation of the stem-cell hypothesis [29] and for exploring FDG-PET (fludeoxyglucose positron emission tomography) dynamics [30] for hypoxia in solid tumours. Crucially, the non-homogeneous oxygen distributions in tumour spheroids have been well studied [31–33]. Research to date shows that cellular oxygen consumption rate (OCR) has a known influence on the oxygen concentration throughout a spheroid, and directly influences the extent of central anoxia and the viable rim thickness by known mathematical relationships [32]; thus measuring these aspects allows an experimenter to determine OCR with relative ease compared with other methods [28].

Such methods and the underlying theory are exceptionally important to understanding the factors that influence tumour oxygen distribution, yet all these methods rely on an implicit assumption of perfect sphericity. There is a clear rationale behind this, as symmetry considerations simplify the problem greatly. Yet in experimental conditions, imperfect spheroids are common, frequently growing as extended ellipsoids. When the eccentricity of these shapes is extreme, an experimenter may reasonably choose to discard them from analysis. But this prompts a question: quite how extreme do such deformations have to be before a spherical assumption breaks down? Presumably small departures from sphericity should not impact analysis, whereas highly eccentric ellipsoids could reasonably be presumed to violate the underlying theoretical assumptions. The question of how such eccentricities might skew analysis of spheroids, and how trustworthy results of such analysis might be, has not yet been considered in the literature, despite its obvious practical importance.

These questions are as of yet unanswered, and are of paramount importance given the growing adoption of spheroids for cancer research, and their utility in estimating OCR [31,32]. Knowing the acceptable limits of eccentricity for spheroid analysis would be of considerable benefit to experimenters, providing error estimates and limits of reliability. A full analytic expression for ‘ellipsoidals’ (analogous to the spherical case) would also be of substantial benefit, allowing the analysis of eccentric shapes and potentially increasing experimental throughput. In this work, we seek to address these issues by deriving an expression for oxygen diffusion in both prolate and oblate geometries. This is contrasted to the spherical case to determine the limits of validity for experimental data, and the implications of this are discussed. A schematic of this is depicted in figure 1.

### 1.1. Spheroids and ellipsoids

The general equation of an ellipsoid is given by
1.1where *a*, *b* and *c* are the major axes' lengths. For an ellipsoid with azimuthal symmetry, *a* = *b*. When all axes are equal (*a* = *b* = *c*), the result is a perfect sphere. To date, this is the only case which has been well studied from a theoretical standpoint [31–33]. While the mathematical treatments to date have assumed spheroids are perfect spheres, the nomenclature ‘spheroid’ still applies to the more general case, including prolate and oblate spheroids. In this work, we broaden the mathematical framework to be applicable to ellipsoids without full spherical symmetry, which are namely

(1)

*Prolate spheroids*. In the case where*c*>*a*, the resulting ellipsoid is an ellipse rotated around its major axis, the line joining its foci. This yields a rugby ball-type shape.(2)

*Oblate spheroids*. Where*a*>*c*, an oblate spheroid results, equivalent to an ellipse rotated around its minor axes. The resulting shape is discus-like.

Examples of these ellipsoids are shown in figure 2. To avoid confusion, we use the term spheroid to refer to an ellipsoidal collection of cells, although we do not limit this to perfect spheres, qualifying with terms ‘prolate’ or ‘oblate’ as appropriate. We use the term ellipsoid to refer to surfaces of iso-concentration of oxygen.

## 2. Model derivation

The full mathematical derivation for oxygen partial pressure in prolate and oblate spheroids is rather involved, and here we shall confine ourselves to stating results with a cursory outline of how they are derived. A full mathematical outline is provided in the electronic supplementary material, appendix, S1. Essentially, we are concerned with solving a steady-state reaction–diffusion problem for oxygen field *P* of the form
2.1where *D* is the oxygen diffusion constant in water (typically *D* = 2 × 10^{−9} m^{2} s^{−1}) and *aΩ* is the oxygen consumption rate in mmHg s^{−1}. This must be solved subject to two crucial boundary conditions, namely that the surface flux and oxygen partial pressure at the anoxic boundary must both be zero. For simple geometries such as perfectly spherical spheroids and cylindrical vessels, symmetry can be exploited to readily yield analytical solutions [33]. In elliptical geometry, the problem is more involved but the basic premise remains the same, and is outlined below. The geometry of the problem is illustrated in figure 3.

### 2.1. Prolate spheroids

In a prolate spherical geometry, we employ the prolate spherical coordinate system, using a geometrically intuitive definition where curves of constant *σ* are prolate spheroids, while curves of constant *τ* correspond to hyperboloids of revolution [34]. This is outlined in detail in the electronic supplementary material, text S1. This yields an analytical solution, which can be converted directly into spherical coordinates to yield
2.2where *f* = *er*_{n} is the distance from the ellipse centre to the foci.

### 2.2. Oblate spheroids

In oblate spherical geometry, a similar geometrical definition exists [35] and can be solved through similar methods, also outlined in the electronic supplementary material, S1. The full solution in spherical coordinates is
2.3Both the prolate and oblate form can be alternatively cast in terms of *r*_{o}, the outer semi-major axis length, if preferable. These forms are also given in the electronic supplementary material, S1.

### 2.3. Ellipsoidal confocality

Analogous to the perfect spherical case, confocal elliptical surfaces in a spheroid are at the same oxygen partial pressure. For confocal ellipsoidal shells, focal length is constant, related to the eccentricity *e*_{c} and semi-major axis of the shell *r*_{c} by
2.4It follows that the innermost (anoxic) and outermost ellipsoidal shells are confocal, thus for a true spheroid *e*_{o}*r*_{o} = *e*_{n}*r*_{n}. Within the bounds of acceptable experimental error, this relationship can be used to determine whether a given spheroid displaying apparent eccentricity is ellipsoidal or not. This is important from an experimental perspective, as sectioning can introduce serious distortions in fixed spheroid sections, or can miss the central axis of the spheroid [32]. In these cases, an ostensible ellipsoidal shape might be observed, but may in fact be a sectioning distortion or off-centre cut. Testing for confocality thus determines the underlying reality.

### 2.4. OCR estimation in spheroids

In the perfectly spherical case, OCR (in mmHg s^{−1}) is related to the anxoic radius *r*_{n} and outer radius *r*_{o} [28,32] by
2.5For a prolate tumour spheroid, it is possible to estimate OCR in a manner analogous to the spherical case by re-arranging the equations for *P*_{P} to yield
2.6Similarly for oblate spheroids, OCR is given by re-arrangement of *P*_{O} to arrive at
2.7One complication that may arise is that it may be impossible to ascertain whether an ellipsoidal spheroid is prolate or oblate. In that case, one can produce a ‘combined’ expression for average OCR by taking the average of equations (2.2) and (2.3) and re-arranging to arrive at
2.8An example of the implementation of these forms including error analysis is included in the code given in electronic supplementary material, S2.

### 2.5. Spherical error metrics

It is worthwhile introducing metrics to quantify how divergent the estimated oxygen profile in a given ellipsoidal spheroid is from a related perfectly symmetric spheroid. A perfect spheroid has radial symmetry, and thus *P*(*r*_{o}) = *p*_{o} at all points. Consider related prolate and oblate spheroids with both semi-major axis *r*_{o} and eccentricity *e*, nested inside a sphere of radius *r*_{o}. We can define the root mean square error (RMSE) by contrasting the expected outer shell partial pressure *p*_{o} with what would be measured for spheroids at *P*_{P}(*r*_{o}, *θ*) and *P*_{O}(*r*_{o}, *θ*), respectively. The RMSE for prolate and oblate spheroids, respectively, is
2.9and
2.10These equations can readily be solved by numerical integration methods, and solutions are demonstrated in the electronic supplementary material, code S2. Percentage error is simply 100(RMSE/*p*_{o}), and thus the variation in RMSE with eccentricity can be readily calculated. The other instance when deviation from spherical assumptions must be quantified is in OCR calculation; for example, when the OCR in a spheroid is calculated assuming the perfectly spherical form in equation (2.5) rather than a more appropriate prolate or oblate form. This might occur when eccentricity is low and the spheroid appears to be entirely symmetric to a first approximation. The distance from the centroid of a spheroid with semi-major axis *r* and focal length *f* to a point on the spheroid at an angle *θ* is given by , and thus the average values for outer radius *r*_{o} and anoxic radius *r*_{n} are given by integrating this over a full revolution, yielding
2.11and
2.12where *E*_{m} is the complete elliptic integral of the second kind. Analogous to the discrete standard deviation, the distance function can be integrated over a full rotation () to yield an expression for standard deviation across these spheroids of
2.13and
2.14Applying the form in equation (2.5) for spherical OCR yields an approximation (which is incorrect when *e* > 0) of
2.15The uncertainty calculation associated with this can be calculated with the variance formula, in this case given by
2.16This is analytically tractable, and is given in the electronic supplementary material, S1. An implementation in several code languages is also provided in the electronic supplementary material, S2.

## 3. Methods

### 3.1. Spheroid oxygen profiles

The models derived in this work were used to create oxygen profiles for spheroids of both prolate and oblate classes, which were contrasted to conventional perfectly spherical profiles.

### 3.2. Quantifying differences between prolate and oblate cases

Prolate and oblate forms have some mathematical differences, as can be seen by inspection of equations (2.2) and (2.3). Whether this difference is experimentally significant is an important question; from a single projection of a tumour spheroid it might be impossible to ascertain whether an experimenter is dealing with a prolate or oblate case. As an experimentalist might not be able to determine whether a given spheroid is prolate or oblate from a single section, quantifying differences in measured OCR under each assumption is an important goal of this work. This was simulated by producing prolate and oblate spheroids with properties as outlined in table 1, and observing the differences in their profiles. In addition, OCR estimates under the ‘wrong’ assumptions were also calculated and inspected. Specifically, the ‘wrong’ assumption occurs when one applies oblate equations for a prolate spheroid or vice versa. OCR was also calculated with the combined assumption (equation (2.8)), which can also be used when the underlying form is unknown. These spheroids were produced with physical properties as per table 1, with OCR given by equations (2.5)–(2.8).

### 3.3. Comparisons with the spherical case

As OCR from spheroids is estimated assuming spherical symmetry, a major aspect of this work was quantifying precisely how close to perfect sphericity spheroids must be so that such an assumption holds, and how departures from sphericity impact estimates of OCR and oxygen profiles. To study this, spheroids with known OCR and varying eccentricity were simulated, and analysed with equations (2.9)–(2.16).

### 3.4. Experimental proof of concept

To date, non-spherical tumour spheroids have been somewhat neglected, and have been frequently discarded from analysis due to their inherent uncertainty. It is thus difficult to find non-spherical tumour spheroid data. A potential example was taken from a previously analysed set of sectioned DLD-1 tumour spheroids [32], dual-stained with proliferation marker Ki-67 and EF5. Spheroids from this set were experimentally determined to have an OCR of 22.10 ± 4.24 mmHg s^{−1}. The sample spheroid was excluded from prior analysis because of its high eccentricity (external eccentricity *e* ≈ 0.66). This was then analysed using methods outlined in this work as a proof of concept to determine OCR, contrasting it with known values and spherical estimates.

## 4. Results

### 4.1. Oxygen distributions in spheroids

Eccentric spheroids were simulated with properties shown in table 1. Unlike the spherical case, oxygen profiles here are not radially symmetric, so profiles were plotted along both the semi-major and semi-minor axis for clarity. Examples of these profiles are depicted in figure 4. Oxygen gradients through spheroids are simulated in figure 5, for both prolate and oblate cases with increasing eccentricity.

### 4.2. Quantification of differences in prolate and oblate spheroids

From two-dimensional sectioning or imagining alone, it can be experimentally difficult to ascertain whether a given spheroid is either prolate or oblate. It is thus important to quantify differences between the ellipsoids. Figures 4 and 5 suggest that prolate and oblate spheroids have broadly similar oxygen profiles until eccentricity approaches unity. Figure 6 depicts internal anoxic radii with eccentricity. These differ only slightly, typically less than or equal to 1% for both major and minor axes for 0 < *e* < 0.9. More interesting perhaps is the variation in OCR estimate with eccentricity under the ‘wrong’ assumption (namely assuming oblate form when the actual entity is prolate or vice versa) shown in figure 7. This suggests strongly that OCR estimates arrived at under the ‘wrong’ assumption are still accurate up until high eccentricity. Even at very high eccentricity, an average value of both incorrect OCRs was extremely close to true OCR. This suggests that incorrectly specifying the type of spheroid should not greatly impact OCR estimates. For improved accuracy, employing the combined estimate in equation (2.8) yields even smaller errors, suggesting this could readily be used by experimentalists without introducing major error even when the underlying form is unknown.

### 4.3. Comparisons with spherical case

Table 2 depicts the impact of assuming perfect sphericity on derived OCR estimates for properties in table 1. At low *e* (typically *e* ≤ 0.3), treating spheroids as perfect spheres yields acceptable accuracy for both OCR estimates. However, as , the reliability of OCR estimates rapidly breaks down, and errors become increasingly large and unreliable.

### 4.4. Experimental proof of concept

For this work, a simple image analysis algorithm was written for the spheroid image, which found the ellipsoid centre and cast best fit ellipses from this position. The analysis algorithm was broadly similar to previously described methods [32], yielding estimates of *e* ≈ 0.66, *r*_{o} = 488.5 μm and *r*_{n} = 400.15 μm. This analysis suggested the best-fit inner and outer ellipses were approximately confocal to within experimental error, as required by the confocality condition in equation (2.4) (*e*_{o}*r*_{o} ≈ *e*_{n}*r*_{n} to within an error of 1.13%). Uncertainty on the lengths of *r*_{o} *r*_{n} were taken from this to be 5.53 μm and 4.52 μm, respectively. Results of this analysis are depicted in figure 8 and table 3, and suggest values in agreement with those previously measured when considered as a prolate/oblate spheroid, and unrealistic values if presumed perfectly spherical.

## 5. Discussion

Analysis of spheroids to date tends to pivot on the presumption of sphericity, as symmetry arguments reduce the complexity required. However, real spheroids tend to depart from perfect sphericity to varying extents. In this work, we provide a metric for determining how reliable the simpler spherical assumption will be as eccentricity increases, as outlined in table 2. The methods outlined in this work can be employed to generate reliable estimates of OCR and oxygen distribution. Another major benefit of this work is that it allows an experimenter to determine OCR even in non-spherical cases when high eccentricity might otherwise render the spheroids in question ill-suited for analysis. Provided the inner and outer ellipses are suitably confocal, the analysis outlined in this work can be employed, and thus should help increase experimental throughput.

As depicted in table 2, presuming sphericity with eccentric sections yields acceptable accuracy when eccentricity is small, but rapidly begins to produce completely unrealistic results for OCR and massive uncertainty. When analysed as either prolate or oblate spheroids, however, OCR estimates are in agreement with the previously measured values. This suggests strongly that, for eccentricity greater than approximately *e* = 0.3, spherical assumptions for OCR cease to be appropriate and ellipsoidal analysis must be employed. While data for this are currently sparse, we were able to demonstrate the principle on the highly eccentric spheroid illustrated in figure 8, determined by image analysis to have *f* = 318.89 ± 2.57 μm. When analysed as a spheroid, OCR estimates were completely unrealistic with huge uncertainty, as seen in table 3. However, when considered as either a prolate or oblate spheroid, OCR measurements were within previously measured values. This is promising, but, as only a single datum point is available, this should be interpreted solely as a proof of concept.

It is worth noting that, from single section images or microscopy, there is no obvious way to ascertain whether a spheroid is prolate or oblate. As the mathematical forms for these are slightly different, this adds an extra uncertainty and prompts the question about which form is preferable to employ. The analysis in this work (figures 6 and 7) indicates that, even if one incorrectly assumes the wrong form, OCR estimates are still very good, with only minimal errors introduced. This holds with only negligible errors until very high eccentricity. The combined OCR form in equation (2.8) yields only negligible error even when the underlying form is unknown. Thus, an experimenter should opt to use this form for OCR estimation when they have no other information on whether the specimen is prolate or oblate, as this does not introduce large errors even at high eccentricity.

While modelling of elliptical oxygen diffusion has the potential to greatly extend experimental throughput, there are a number of scenarios in which an ostensible eccentric spheroid might not be what it appears. For fixed and sectioned spheroids, the act of sectioning itself can be enough to induce substantial deformations, stretching it along a particular axis. Ostensibly, the resultant shape might appear ellipsoidal, but is in reality a warped spheroid, and cannot be reliably analysed with the methods outlined. Such an example is show in figure 9, for a spherical spheroid sheared along an axis. From the mathematics established in this work, we can distinguish between true spheroids and warped spheroids—if the inner and outer ellipses are not confocal (*e*_{o}*r*_{o}≠*e*_{n}*r*_{n}) then the shape is a warped spheroid, and should be discounted from analysis, as per figure 1. Crucially, figure 9 demonstrates that warped spheroids can only satisfy the ellipsoidal confocality condition under two circumstances: either when its eccentricity is 0 (a perfect sphere), or the non-physical situation when *r*_{o} = *r*_{n}. Thus, a sheared spheroid in one direction will never come close to satisfying the ellipsoidal confocality condition. In practice, all experimental work comes with inherent uncertainty, so *e*_{o}*r*_{o} ≈ *e*_{n}*r*_{n} within the bounds of image analysis uncertainty is sufficient to determine whether a spheroid can be treated as an ellipsoidal case.

There is a more subtle issue with sectioned spheroids, which becomes even more crucial with sectioned eccentric spheroids. Analysis relies on a section through the central axis of the spheroid. In the perfectly spherical case, if the section is off-centre, the net result will be two concentric circles but with a misleading ratio, rendering any OCR calculation derived from this suspect. By contrast, any plane through an ellipsoid produces an ellipse, but if these cuts are off the central axis then the inner and outer ellipses will no longer have a common centre, and will not be confocal. In this regard, determining an off-centre ellipsoid section is relatively straight-forward. A proof of this is provided in the electronic supplementary material, S1.

The theoretical analysis outlined here presents biological investigators with new methods for extending spheroid analysis, and the means to interpret data which depart from sphericity. It also establishes uncertainty bounds on existing spherical analysis techniques, and methods for determining OCR and oxygen distribution in tumour ellipsoids. The findings of this work will increase experimental confidence with tumour spheroids, and have the potential to substantially increase experimental throughput, improving our insights on everything from tumour hypoxia to drug delivery. Such an approach is imperative if we are to fully exploit this unique experimental tool, and ultimately marshal the new insights obtained towards better cancer treatment and diagnostics.

## Data accessibility

All data are available with the paper. Mathematical derivations and related theorems are contained in the electronic supplementary material, S1. Sample implementations in Matlab, Octave, Mathematica and Excel spreadsheet are included in the electronic supplementary material, S2.

## Authors' contributions

D.R.G. and F.J.C. conceived the theory outlined in this work; D.R.G. derived the model, performed the analysis and wrote the manuscript. Both authors reviewed the manuscript and authored supplementary texts.

## Competing interests

We declare we have no competing interests.

## Funding

The authors acknowledge Queen's University Belfast for funding the CAIRR initiative. D.R.G. also thanks Cancer Research UK for its support.

## Footnotes

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.4183505.

- Received April 16, 2018.
- Accepted July 23, 2018.

- © 2018 The Authors.

Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.