Intravital mesoscopic fluorescence molecular tomography allows non-invasive in vivo monitoring and quantification of breast cancer growth dynamics

Preclinical breast tumor models are an invaluable tool to systematically study tumor progression and treatment response, yet methods to non-invasively monitor the involved molecular and mechanistic properties under physiologically relevant conditions are limited. Here we present an intravital mesoscopic fluorescence molecular tomography (henceforth IFT) approach that is capable of tracking fluorescently labeled tumor cells in a quantitative manner inside the mammary gland of living mice. Our mesoscopic approach is entirely non-invasive and thus permits prolonged observational periods of several months. The relatively high sensitivity and spatial resolution further enable inferring the overall number of oncogene-expressing tumor cells as well as their tumor volume over the entire cycle from early tumor growth to residual disease following the treatment phase. Our IFT approach is a promising method for studying tumor growth dynamics in a quantitative and longitudinal fashion in-vivo.

T he current frontiers in the battle against breast cancer are shifted towards early detection on the one hand and understanding mechanisms that lead to recurrent tumors on the other hand. While therapeutic improvements in combination with early tumor detection enabled through the introduction of mass mammographic screening programs has led to a substantial increase in breast cancer survivors over the last few decades 1 , a detailed, longitudinal insight into early tumor initiation events under physiological conditions is still lacking 2 . A further major problem in modern breast cancer care is constituted by the increase of breast cancer related death due to recurrent disease 3 . These incurable recurrences are caused by treatment refractory cells that stay dormant over years 4,5 . Similar to the situation at the tumor induction phase, these small amounts of residual cells are hard to trace and molecular markers of this population need to be established to understand potential points of therapeutic interference 6 . Such information may be difficult to decipher solely by studying cancer in the human population, due to the impossibility to follow scarce cells in these cancer stages for mechanistic analyses in an enormously heterogeneous patient population.
To this end, preclinical mouse models can provide valuable resources for cancer prevention and treatment strategies, since they enable both longitudinal and molecular analyses of precancerous and cancerous lesions in defined genetic contexts and in environmentally controlled conditions, which would be impossible to accomplish by studying human cancer 7,8 . Preclinical breast cancer mouse models, such as transplantation based mouse models or genetically engineered mouse models (GEMMs), have been successfully used to investigate tumor biology via a range of imaging approaches 9 and therapeutic strategies 10 . Moreover, GEMMs based on conditional gene targeting wherein genes of interest are inactivated (or activated) in a temporal and tissuespecific manner offer advanced tools to study cancer with the inclusion of reporter alleles to perform in vivo imaging 11 . Finally, recent advances in mammary tumor biology with respect to transplantation 12 and intraductal injection techniques 13 permit the study of cells derived from above-mentioned models as well as genetically engineered lines (Fig. 1a, b).
This extensive collection of available preclinical breast cancer models represents a powerful toolset to explore mechanisms on tumor evolution and therapy success. Despite their unique advantages, current studies employing preclinical mouse models do not harness their full potential. This is because currently employed methods to visualize and assess tumor progression have important limitations with regards to either sensitivity, resolution, imaging field-of-view (FOV), and depth or the ability to monitor disease over prolonged timeframes. In particular, histological approaches, while providing valuable quantitative, high-resolution information on tumor morphology and abundance, are end-point methods and thus prohibit reconstruction of individual, longitudinal tumor development. This ability, though, is critical in assessing tumor populations during early growth phases as well as during treatment to understand dynamics of shrinkage and evolution of resistant disease. Intravital, longitudinal imaging techniques, on the other hand, can give valuable information on the temporal evolution and dynamics of cell populations 14 . But because of the thickness 15 (~300 µm) and highly scattering nature of the mouse skin and tumor tissues, optical access for standard light microscopy methods such as confocal or two-photon has to be guaranteed by either surgical implantation of so-called imaging windows 16,17 or by applying windows on surgically exposed organs [18][19][20][21][22] . In addition to the limitations in imaging FOV and depth, these surgical procedures, however, typically trigger inflammatory responses that may alter the tumor microenvironment 16 , are limited in their respective size and bear major burden for the animal when utilized over extended period of time. Furthermore, as these windows are static, they are not fully compatible with spontaneous mouse breast tumor models, as these cases require dynamic and flexible repositioning of the imaging field-of-view (FOV) in order to monitor stochastic tumor emergence over 10 mammary glands.
An alternative method for non-invasive monitoring of tumor burden and screening for drug efficacy in cancer mouse models is bioluminescence imaging (BLI). Although its non-invasive character and high sensitivity are promising features 23 , the spatial resolution achievable in 3D reconstructions is typically very poor (~mm's) 24 . Furthermore, the bioluminescent signal mainly depends on two components: the substrate/enzyme amounts and their interaction, as well as the accessibility of oxygen in the microenvironment. Both parameters are variable and challenging to assess independently, therefore the normalization of bioluminescent signal required for extracting quantitative information has remained challenging. While quantification attempts were shown in 2D imaging 25 , tomographic studies are predominantly qualitative 26 to date, and have instead capitalized on the highthroughput nature and sensitivity of the method. Furthermore, BLI requires injections for every imaging experiment in a longitudinal time-series, and thus puts additional burdens on the animal.
Besides the above-mentioned fluorescentand luminescent-based techniques, other mature technologies are also being used for monitoring tumor development. These conventional techniques include computed tomography (CT), magnetic resonance imaging (MRI), and ultrasound (US) and typically only deliver unspecific structural information with limited spatial resolution 27 . Despite recent advancement that enable cellular specificity through the use of engineered molecular constructs 28 , the achievable spatial resolution is again limited to a few hundred microns. Furthermore, for CT and MRI, achieving high spatial resolution requires prolonged imaging sessions that may last over several hours, which puts additional burden on experimental animals. Optical spectroscopy techniques can be regarded as an alternative non-invasive, fast, and low-cost approach 29 . As these methods quantify the spatiotemporal changes of physiological parameters (e.g. oxy-/deoxy-hemoglobin concentration, oxygen saturation etc.) of the tumor environment 30 , they however only provide surrogate (e.g. metabolic) information about the state of tumor cells. Another recently emerging technique is photoacoustic imaging (PAI) which is capable of detecting absorbing contrast agents deep (~1 cm) in scattering tissue with mesoscopic resolution (a few 100 µm). PAI is inherently sensitive to chromophore contrast so it is suitable for imaging vascularization and their oxygenation state via intrinsic spectral absorbance differences of (oxy-/deoxy-) hemoglobin. Moreover, the recent introduction of gene reporters for molecular imaging in Photoacoustic Tomography (PAT) offers the possibility of non-invasive imaging with cellular specificity and high sensitivity (~100-1000 cells) 31,32 . Although multispectral PAT is rapidly being developed, current PAT implementations involve costly and custom-made components and its operation entails complicated post-processing, thus limiting its use to a handful of labs worldwide.
Here, we present an experimental approach that aims to address many of the above-mentioned shortcomings. Our platform, named IFT, utilizes an optical tomography approach known as mesoscopic fluorescence molecular tomography (MFMT) [33][34][35][36][37][38][39] . More recently, MFMT methods have been applied e.g. to image the distribution of dyes in mouse brains ex vivo 40 and in vivo 41 , and cellular distributions in engineered tissues 42 or in vitro models 39 . However, applications to tumor imaging have so far been restricted to single timepoints 36 and no long-term longitudinal observations or quantifications of cell numbers have been performed. In our work, we include a number of key improvements that enable quantitative and longitudinal studies of . c Principle of IFT method: (i) focused laser beam excites fluorescence inside the sample volume by raster scanning while a detector array simultaneously images the scanned area. For each scan location, multiple-scattered emission light is collected by the detector array (sCMOS camera). (ii) After a full scan of the sample, each detector yields an image, with all detector images together constituting the raw data (measurement). (iii) A GPU-accelerated Monte-Carlo algorithm simulates the optical propagation of scattered light inside the sample volume. (iv) The raw data and model matrix are fed into a reconstruction algorithm, to extract the 3D distribution of the fluorescence emission. d A CAD representation of the IFT imaging setup. An excitation laser (e) enters the system through a standard filter cube (FC). The laser is raster scanned by a 2D galvanometric mirror (GM) and imaged onto the sample by a scanning lens (SL). The emitted fluorescence light is de-scanned, separated from the excitation light by the FC, and imaged onto the camera (C) by an imaging lens (IL). For in vivo experiments we use an anesthesia system (AS) for narcosis and monitoring of physiological parameters while the region of interest is stabilized with a custom, non-invasive imaging window (NIVOWzoom-in). e Mouse skin position as a function of time as measured by OCT, showing tissue motion along z and x axes without (left) and with (right) NIVOW.
in vivo tumor development in inducible mouse models. First, inspired by Looney et al. 19 , we introduce a modified, non-invasive imaging window that minimizes motion artifacts induced by breathing of the animal and ensures repeated imaging of the same location throughout the study with high stability. Second, subtle changes to the opto-electronic design of our platform compared to previous implementations 39,43 enables larger scan areas (FOVs) at twice the data acquisition rates of earlier versions. Third, after careful characterizing our IFT tomographic resolution and sensitivity, we calibrate its signal with both in vitro phantoms and histological tumor samples, thereby establishing a linear relationship from which we can infer the number of fluorescently labeled, i.e. tumorigenic, cells with high (~200 µm) mesoscopic resolution in 3D. Previous works have monitored tumor xenografts at single timepoints using MFMT 36,[43][44][45] , however, the MFMT signal was never before put into the context of H&E tumor tissue analysis to calibrate the output in quantitative terms. Applying this workflow, here we report -for the first time to the best of our knowledge -longitudinal MFMT imaging of the entire tumor/treatment life-cycle spanning several months in vivo.
We show the capability of IFT to monitor tumor progression over several months in a non-invasive, longitudinal fashion by applying it to a preclinical mouse model of tractable breast cancer. In particular, we track cell number as well as tumor volume and geometry through an entire cycle of tumor growth and tumor regression upon silencing of the oncogenic signal for early (nonpalpable stages) thus relatively small tumors. Here, the high sensitivity of our technique enables tumor identification long before they become palpable, thus making it possible to investigate aspects of early tumor growth and also treatment effects monitoring survival of refractory subpopulations, while the high 3D spatial resolution enables assessing subtle changes in tumor geometry over time.

Results
Principles of intravital fluorescence tomography. Intravital Mesoscopic Fluorescence Molecular Tomography (i.e. intravital MFMT, henceforth IFT) excites fluorescence via one-photon absorption and collects multiple-scattered photons on a 2D detector array (Fig. 1c). The data collection principle is conceptually similar to previous implementations of MFMT 39,46 . An excitation laser beam (561 nm) is focused on the sample and raster-scans an excitation spot across the sample surface. For each scan point, a camera records a wide-field image of the sample, containing all the backscattered/fluorescent photons. Each camera pixel (detector) collects photons at a varying distance away from the laser (source; Fig. 1c, i). With increasing detector distance, the collected photons predominantly originate from increasing depths in the sample, thus the resulting camera image carries intrinsic depth (axial) information. In order to reconstruct a meaningful 3D image from the raw 2D data sets, IFT is based on a so-called forward model (Jacobian matrix), which is calculated by simulation of photon propagation inside the scattering tissue of interest via a GPU-accelerated Monte-Carlo method 38,47 ( Fig. 1c(iii) and Methods section). The optical properties of mammary gland tissue and skin were measured by a home-built dual integrating-sphere setup 48 (Supplementary Note 1 and Supplementary Fig. 1 and 2), and the so obtained values were used as a starting point for the calculation of the Jacobian matrix, which was further optimized to ensured accurate and reliable 3D reconstructions. The Jacobian was then used to retrieve the 3D distribution of the fluorescence emitters from the experimental imaging dataset by running an Lp-norm reconstruction algorithm 49,50 (Fig. 1c(iv) and Methods section).
Non-invasive imaging window. Intravital imaging approaches typically involve the surgical implantation of an imaging window which can trigger inflammatory responses and bear major burden for the animal when utilized over extended time-frames [16][17][18][19][20][21][22] . Moreover, one recurring challenge in intravital imaging is breathing-induced tissue motion, which commonly leads to severe artefacts when not properly accounted for. While the image acquisition can in principle be synchronized or gated to the periodic movement of the body cavity 51 , such techniques are cumbersome and require specialized equipment, are prone to error and typically result in prolonged imaging times as the tissue is at rest only during a small fraction of the breathing cycles. To overcome these challenges, we designed a non-invasive, vacuum-operated stabilization window (NIVOW), inspired by Looney et al. 19 . It comprises a round metal ring (22 mm inner diameter) sealed with a coverslip that can be placed on the mouse skin at the center of the imaging FOV (Fig. 1d). Vacuum can be applied through two small openings, which firmly attaches the skin to the coverslip, while the window housing is secured to the optical table via suitable arms. This approach allowed us to leave the mouse skin intact (i.e. no surgical intervention) and enabled repeated and flexible positioning for imaging. Since the tumor mass is attached to the skin and both, tumor and skin, can be lifted away from the body wall by applying minimal suction, we do not expect this procedure to cause recognizable tumor shape deformation. Furthermore, the vacuum stabilizes the region of interest and its underlying tissue substantially, both laterally as well as axially (Fig. 1e). Quantitatively, the axial motion of skin and subcutaneous tissue layers decreased from over >100 µm to~1 µm peak-to-peak, as measured at high speed by optical coherence tomography (OCT; Supplementary Note 2 and Supplementary Movie 1 and 2). This residual motion is negligible compared to our spatial resolution and the coverslip further ensured a 'flat boundary' condition as required by the Jacobian model matrix. Thus, our NIVOW is an instrumental part of the IFT system for achieving high-fidelity reconstruction. Furthermore, it enables flexible repositioning of the imaging location. As the site of tumor onset can shift in the abdominal mammary gland area by several millimeters, having an adjustable region of interest brings a critical advantage over traditionally used, fixated imaging windows 16 . Most importantly, however, no surgical intervention is required for placing the window and thus (i) no inflammatory signal can interfere with the tumor biology under investigation, and (ii) monitoring can be performed over the entire life-time of the animal (see Supplementary Movie 3).
IFT characterization (resolution and sensitivity). Our aim is to capture the process of breast cancerous tumors in vivo in its entirety: from the first formation of a solid tumor mass by a cluster of cells to the stage where the tumor becomes barely palpable (typically >5 mm) until the stage of very small populations of residual tumor cells following therapy. For this our IFT imaging system must display the required spatial resolution as well as sensitivity. Initial tumor populations are commonly assumed to comprise~1 × 10 4 cells, which hence marks our sensitivity requirement. Considering the average diameter of a mammary gland tumor cell (10-15 µm), these smallest units of tumor populations approximate to 0.01 mm 3 volume, hence necessitating a similar spatial resolution on the order of~250 µm. In order to quantify the spatial resolution of IFT, we performed a terminal experiment and placed individual, 50 µm fluorescent polymer microspheres (Cospheric, USA) under the mouse skin at a depth of~750 µm (as measured with a caliper). We then imaged these beads with our IFT system (Fig. 2a). Reconstructing the acquired IFT image dataset indeed yielded single, bead like objects, which we used to characterize the point-spread function (PSF) of our system (Fig. 2b). In contrast to previous studies which often characterized spatial resolution by reconstructing a single focal object, here we purposely chose to reconstruct two objects in close proximity. Our characterization yielded a lateral and axial resolution of~200-220 µm, and~250 µm, respectively, thus sufficient to resolve initial tumor clusters of~10,000 cells. While resolution and sensitivity are generally expected to decrease with depth in fluorescence tomography 34,52 , the changes over our imaging and reconstruction depths are within the limit of our discretization (i.e. 1 voxel, laterally) and thus assumed to be negligible (also see Methods section).
Reconstruction algorithm linearity. One of the main aims of our study which goes beyond the capability of other established cancer imaging methods is to quantitatively measure the number of oncogene-bearing cells that contribute to the signal in IFT. A critical prerequisite for this is to show that the intensities of the IFT reconstructions scale linearly with the underlying fluorescence intensity. To investigate this relationship in vitro, we engineered a scattering tissue phantom (Intralipid 10%) comprising fluorescein filled tubes (~3 µL) with three different, physiologically relevant, concentrations (25, 2.5, and 0.25 µM) at a depth of 1.3 mm (Fig. 2c). These depth and concentration range was chosen since it represents the typical tumor location and fluorescence signal ranges. We then collect the data for these concentrations separately and normalized signals based on excitation power. IFT reconstructions (Fig. 2d) of these phantom experiments delivered signal levels that scaled linearly (R = 0.996) with actual fluorescein concentrations (Fig. 2e). We did not assess linearity outside this concentration range, as we could not reliably detect lower concentrations and never encountered higher signal levels in vivo, respectively.
Calibration of IFT cell count. Since in vitro phantoms cannot replicate the real, in vivo tissue heterogeneity, we next proceeded to calibrate our IFT signal values with absolute cell quantities obtained through tumor histology (Fig. 3). The presence of strong heterogeneity in nuclear H&E staining made it difficult to count nuclei with a standard cell counting software, therefore we devised a custom image processing pipeline (Fig. 3a,b). Here, we first restricted our analysis to a small sub-area of the tumor site that displayed homogenous H&E staining, and manually identified tumor clusters for automatic cell counting. This yielded crucial tumor cell density information (#cells/area). In parallel, the entire H&E slice was thresholded to separate tumor cells from other, stromal areas. Then, by measuring the respective entire tumor area and utilizing the cell density information, the overall number of cells pertaining to a given tumor section could be calculated. Likewise, this calculation was repeated for all histological sections, yielding an overall cell number for a 3D tumor. This pipeline was  Figure 3c shows histological (ground truth) number of cells plotted against the integrated IFT reconstructed signal for each tumor. We found the resulting graph to display a linear relationship with a fitting performance of R 2 = 0.94 and RMSE = 4.78, thus demonstrating the ability of our IFT system to faithfully capture the amount of fluorescent, i.e. oncogene bearing, tumor cells from the reconstructions. The fit to this data constituted our IFT calibration and was hence used in all experiments to convert reconstructed IFT signal to actual cell quantities. Furthermore, we note that the smallest tumor investigated during this calibration study was histologically measured to only contain~10 5 cells. As we were also able to identify and reconstruct this small tumor mass via IFT this number also represents a conservative estimate of our in vivo IFT sensitivity. In fact, the spatial distribution of tumors cell count (C sub ). Overall tumor cell number was then estimated as indicated. Note that this sub-area analysis was repeated for every slice of each tumor. Scale bar, 100 µm (zoom-in 20 µm). c Integrated IFT reconstruction signal values for (n = 8) independent tumors of varying size and intensity plotted against 'ground-truth' cell count as obtained by H&E analysis. Linear fit to the data constitutes the IFT cell count calibration and was used for subsequent experiments. d Exemplary 2D IFT cross-sections and corresponding H&E slices as well as 3D IFT reconstructions for two representative tumors in C (also see Supplementary Fig. 3 for a slice-by-slice visualization of IFT 3D reconstructions). Scale bars 1 mm (left) and 2 mm (right, 3D reconstruction). (Fig. 3d) indicates a detectable signal levels on the order of a few thousand cells. Finally, IFT reconstructions and histological slices of the same tumor show considerable morphological overlap, further corroborating the capability of IFT to faithfully capture the 3D tumor geometry and volume (Fig. 3c, right). Also see Supplementary Fig. 3 for a slice-by-slice visualization of IFT 3D reconstructions.
Longitudinal molecular imaging of tumor development. After characterizing our IFT imaging system, we applied it to monitor mammary gland tumor development longitudinally in a proof-ofprinciple experiments (Fig. 4). For this, mice (n = 2) with primary tumor cells implanted into their fat pad area were put on Doxycycline (DoX) for 4-5 weeks, throughout which the tumorigenic cells continued to grow (see timepoints t 1 -t 2 in Fig. 4). This suspected growth phase was clearly recapitulated by a strong (~10-fold) increase in IFT signal and hence tumorigenic cell quantity. Likewise, after the removal of DoX from the food diet, the cell numbers rapidly declined within a week (t 3 ), very much similar to the dynamics observed in the original transgenic model 6 , and returned to or below base level within~40-50 days (t 4 ). At this point (t 4 ) the remaining, weak IFT signal can be attributed to scar tissue of the residual disease, which we therefore consider as the baseline (autofluorescence) signal. From these reconstructed data, we can also estimate the sensitivity of our system, which is on the order of 2-3 × 10 3 cells per voxel. These tumor induction and regression dynamics could be replicated in additional animals (n = 4, see Supplementary Fig. 4). Remarkably, although we could observe a clear dependence of cell quantity on DoX, overall tumor volume did not change substantially with the DoX cycle or cell quantity as assessed by IFT.

Discussion
In this work we have presented a proof-of-principle demonstration of a non-invasive and quantitative approach to image mammalian tumor development with molecular contrast at high mesoscopic resolution and sensitivity. Our intravital fluorescence tomography system uses a non-surgical and re-useable imaging window and thus permits prolonged observations of tumor development over several months and with arbitrary frequency, and is thus able to capture in principle the transition from tumor induction to tumor progression all the way to tumor regression upon treatment. One advantage of our method is the capability to resolve relatively small tumor cell masses in 3D, comprising only few hundreds of µm and~10 4 cells overall, which would normally be missed by palpation or other established methods such as ultrasound and MRI 28,53 . From a biological point of view, these capabilities could enable a close evaluation of nascent stages of tumor induction, which could be detected and 3D-reconstructed upon IFT recording. Depending on the 3D tumor geometry changes and the observed growth rate between measurements, tumors can be isolated for culture and analysis at informed, early time points that are normally not accessible for non-palpable tumors. Similarly, the sensitivity of the IFT measurements could aid to examine stages of tumor regression, opening the possibility to capture treatment of refractory tumor clones that are not palpable anymore (e.g. time point t 3 in Fig. 4). These normally inaccessible cell populations, also termed minimal residual disease, are known to be responsible for incurable tumor recurrences 4 . Hence, a handle to follow these cells and more resistant clones during the process of tumor treatment and the chance to isolate and study them in detail in vitro would yield further insights on late stage tumor biology 5,8 .
Furthermore, a particular strength of our IFT technique is the ability to infer quantitative information about the underlying amount of fluorescently labeled (i.e. oncogenic) cells at relatively high accuracy. This capability sets our approach apart from other non-invasive imaging methods such as MRI, ultrasound, CT, BLI, or photoacoustics, which can provide only information about overall tumor volume at mesoscopic resolution.
A critical step in rendering IFT quantitative is the calibration to histological ground-truth data, which requires additional animals and sophisticated, time-consuming analysis. However, once calibrated, the conversion factor is expected to remain reasonably constant for set of animals with similar genotype, gender, and age etc. 15 , and thus can be used across animals and experiments, provided that the scattering properties of the tumor and surrounding tissue and the expression of molecular fluorescence reporters does not change substantially between the calibration and experimental batches. Here, our recent work 54 suggests that fluorophore expression and thus brightness of individual oncogenic cells stays constant over time, and is not affected by the presence of DoX. Thus, an increase in overall IFT fluorescence can indeed be attributed to an increase of oncogenic cells within a given volume. Other experimental uncertainties entails cell counting variabilities due to difference in histological staining as well as spatial inhomogeneity in tumor cell sizes and thus densities between the analyzed and extrapolated histological areas used for calibration (c.f. Fig. 3a, b). We estimate the overall uncertainty in absolute cell numbers to be on the order of 20-30%, which is based on extrapolations of extremal densities found in a representative tumor (see Supplementary Table 1). Furthermore, tumor compression and shape changes cannot be completely ruled out with our vacuum-operated window, however, such effects are not expected to change the cell number nor total volume and thus do not alter our main measurement outcomes.
An asset of IFT is its capability to provide information on the 3D geometry of tumors early in progression and during unpalpable stages of tumor regression will enable assessments of invasive or satellite states of the tumor and can be correlated with in depth histological and functional analysis upon informed resection from the animal (Fig. 3d). Here, one unexpected result of our study was the fact that tumor volume and (oncogenic) cell quantity are not well correlatedthis could be because in different stages of tumor development, the density of similarly sized tumors can vary substantially due to the different amounts of non-fluorescent, stromal cells that infiltrate the tumors. Therefore, also the overall density and thus the number of fluorescence tumorigenic cells within a given volume can vary (see Supplementary Fig. 5). The assessment of quantitative tumor volume in IFT depends critically on thresholding parameters, and can in principle recapitulate true tumor volume as assessed by H&E (see Supplementary Fig. 6). Here, our unbiased thresholding approach based on mutual information metric can provide more accurate volume measurements compared to standard techniques in the field 49 (see Methods section, Supplementary Note 3 and Supplementary Figs. 7-10). However, in absence of such H&E ground truth information, any absolute volume measurement has to be treated with caution, albeit relative changes over time can be reliably observed (see Fig. 4). In cases where information on quantitative tumor volume is desired, one has to extract the exact optical properties of the underlying tissue for reconstruction with using techniques suggested in the literature 55 or can resort to alternative methods such as US or PAT.
Unlike imaging techniques that achieve higher spatial resolution and molecular sensitivity such as confocal or multiphoton microscopy, IFT does not require surgical intervention to gain optical access to the tumor, and is capable of imaging at a much greater tissue or tumor depth. In contrast to BLI approaches, IFT does not necessitate external substrate injection during the imaging process, and thus IFT can collect data with high consistency for months without compromising the animal well-being. Combined with the relative simplicity and low-cost of the experimental apparatus, we expect IFT to become an attractive method for tracking in vivo tumor progression in deep mammalian tissues in a non-invasive and longitudinal fashion. In the future, the fact that IFT can in principle record fluorescence at multiple wavelengths quasi-simultaneously should enable to study the spatial and temporal dynamics of different clonal tumor populations, their interactions as well as their effect on overall tumor burden. In addition, further functional reporters, for instance indicating apoptotic onset at a different wavelength 56 , could be employed to extract additional mechanistic information over the tumorigenesis process.
Taken together, our work and proof-of-principle demonstrations show the strengths of preclinical mouse models when combined with the non-invasiveness, high resolution, and sensitivity of lightbased imaging methods. This might enable future studies on therapeutic interference with critical early stages of tumor development as well as on treatment refractory tumor clones.
Statistics and reproducibility. No prior assumptions were made regarding effect sizes. All relevant and possible data points were used for statistical comparisons when applicable. Biological experiments were performed independently of each other, including media preparation. No acquired data was excluded from the analysis. No randomization was applicable. Mouse strains utilized in this study: RAG1 (−/−), as immunodeficient mouse receiver strain 57 (C57BL/6 background), and the trackable tumor donor strain TetO-cMYC/TetO-Neu/MMTV-rtTA/H2B-mCherry [FVB background] (i.e. tumorigenic cell nuclei are labeled with the fluorescence protein mCherry). In order to generate the donor strain, the tri-transgenic mouse line TetO-cMYC/ TetO-Neu/MMTV-rtTA 6,58,59 was bred with the R26-H2B-mCherry strain 60 .
In this mouse model tissue-specific induction of oncogene expression results in the rapid onset of invasive mammary carcinomas. These tumors display strong oncogene addiction; hence, following oncogene inactivation they completely regress to a non-palpable state. The regulation of the oncogenes is achieved by controlling them with the TeT-on system. In short, the Tet-On expression systems is a binary transgenic system in which expression from a target transgene is dependent on the activity of an inducible transcriptional activator. In the Tet-On system, expression of the transcriptional activator can be regulated both reversibly and quantitatively by exposing the transgenic animals to varying concentrations of doxycycline (Dox).
All efforts were made to minimize the amount of animals used in accordance with Russell and Burch's principle of (3Rs) reduction and highest ethical standard. In total 30 mice were used for the experiments lined out in this manuscript which includes mice for preliminary tests on anesthesia conditions and fat pad injection technology as well as injected cell numbers needed for tumor induction. Animals were kept on a 12-h light/12-h dark cycle, with constant ambient temperature (23 ± 1°C) and humidity (60 ± 8%), supplied with food pellets (for tumor induction, pellets contained Doxycycline hyclate, 625 mg/kg; Envigo Teklad) and water ad libitum.
Cells transplantation. Mammary glands from the donor strain (female virgin mice 7-9 weeks old) were harvested, singularized and seeded on a 2D collagen coated sixwell plate (BD, Cat. # 356400). Cells were fed with mammary gland specific media (Promocell, Cat.#21210), plus Doxycycline hyclate (Sigma, D9891) overnight. The next day, tumor cells were extracted, counted, and resuspended in a solution of PBS and Matrigel (Corning, 356231). A mixture of 150,000 cells in 40 µl of media were injected into the donor strain using the 50 µl syringe (Hamilton, 705 N ga22s/51 mm/ pst2). Female virgin RAG1 animals (from 4 to 8 weeks), were fed with Doxycycline diet (Envigo, TD.01306) at least 2 days prior to transplantation and for as long as tumor growth was desired. Animals were anesthetized using Xylazine/Ketamine before the surgery procedure. Closure of the wound was done by manual suture and no further post operational issues or inflammatory responses were observed.
Optical setup. A schematic overview of the experimental arrangement is given in Intravital fluorescence tomography. In order to reconstruct a 3D image from a set of 2D images, IFT proceeds as follows. First, a so-called forward model (Jacobian matrix) is calculated, in which photon propagation inside the scattering tissue is simulated via a GPU-accelerated Monte-Carlo method 38,47 (MCX, 2019.4, GeForce GTX 1080, 8MB) (see Fig. 1c(iii)). The Jacobian matrix was computed with 10 8 photons and a 0.25 × 0.25 × 0.1-mm 3 voxel size. To mitigate the degradation of resolution along the z-axis with depth, we used a smaller voxel size along this axis. In principle this could also be done for the lateral voxel size, however, then the reconstruction algorithm would be computationally costly and could not be performed on a desktop PC but require more advanced computational infrastructure. On our desktop PC (Intel Core i7-6800K @ 3.4 GHz) computation of the Jacobian matrix took~3 min and total 3D volume reconstruction of the in vivo experiments (10 × 10 × 2 mm 3 volume) took 10 min on average. Note that the Jacobian matrix only needs to be generated only once. Our Matlab code for 3D reconstruction is available on Zenodo 61 .
The integrating sphere experiments (Supplementary Note 1 and Supplementary  Fig. 1 and 2) showed that any differences in optical properties between the excitation (561 nm) and maximum emission wavelength (610 nm) is <5%. Therefore, we used the same optical property values for the Jacobian simulations. For the majority of the results presented here, the following parameters were used: μ s = 20 mm −1 , µ a = 0.1 mm −1 , g = 0.85. We used the same optical properties for skin and tumor tissue in our IFT reconstruction, which is a common practice in the absence of a priori data 36 . For solving the inverse problem of the image reconstruction, the following cost function was used, which delivered an optimized solution for the target volume x: minkA m n x n 1 À b m 1 k 2 2 þ λkx n 1 k 1 ; Here A is the Jacobian matrix and b the experimental imaging dataset. The regularization parameter, λ, is a penalizing term which acts as a filter for x to prevent noise from dominating the reconstruction. An initial value for λ was found by utilizing the L-curve method 62 . Then, a fine tuning was conducted based on the output of mutual information 63 based thresholding (Supplementary Note 3 for details). The regularization parameter (λ) was chosen based on the highest similarity value. This procedure allowed us to account for slight differences in noise levels between individual datasets, while minimizing any user bias in the reconstruction. Before 3D reconstruction, the collected raw images were sub-divided into smaller images depending on the desired source-detector (S-D) separation, which is proportional to the targeted imaging depth and inversely proportional to both lateral and axial resolution. In our work, the chip size of the sCMOS camera allowed us to place detectors up to 3.5 mm away from the source (aligned to the center of the camera). With these parameters, the system is, in principle, capable of reconstructing deeply seated fluorescent emitters from up to 3 mm below the surface, however, in our study we limit our maximum depth to 2 mm. Reconstruction parameters except the regularization parameter (e.g. tolerance and maximum iteration) were kept constant across all experimental time-points. For 3D visualization and assessment of object (tumor) volume while preserving the maximum information on tumor content, we devised the following, non-biased thresholding procedure which is explained in more detail and a step-by-step manner in Supplementary Note 3 and Supplementary Figs. 7-10: First, system based noise in the raw data is removed by taking out pixel values smaller than the median intensity over all 80,688 data points (for 41 × 41 scan positions and 48 detectors). Then, in order to mitigate the diffused signal halo effect, we thresholded the center detector image (i.e zero source-detector separation) by keeping only fluorescence source signal between the maximum intensity and two standard deviation below the mean signal. Note that this is a purely empirical estimation and that this value would need to be adapted in other experiments that utilize different fluorescence reporters or intensities, based on their signal distribution. Next, the raw data was thresholded in an iterative fashion, and compared with the maximum intensity projection (MIP) of the 3D reconstruction in each iteration (see Supplementary Movie 4). To assess overlap, we computed the mutual information (MI) between the raw and the thresholded reconstructed data and chose the threshold with the maximum MI. The MI-based thresholding approach allowed us to minimize user bias when handling data set without a priori data about the size of the target. Moreover, our method allows to also appropriately reconstructing small, scattered signals that belong to tumor protrusions and satellites, which might be missed by a fixed (e.g. 50%) threshold parameter which is common practice in the field of fluorescence tomography 43,64 .
IFT imaging procedure. Mice were imaged according to approved EMBL IACUC protocols. Each animal was anesthetized for the duration of the imaging via Isoflurane. Once heartbeat and breathing stabilized after the initiation of the anesthesia, a topical hair removal cream was applied covering an area of~5 cm 2 revealing the abdominal mammary gland. The imaging procedure consisted of an initial quick scan over the entire area of interest (4th and 5th mammary gland) to identify locations of fluorescence and thus potential tumor initiation. Once identified, we applied the non-invasive imaging window onto the area (Supplementary Movie 3) and acquired the raw IFT data. Typically, this entire process took 15-20 min, of which the raw data acquisition time only took <20 s. During the initial stages of tumor induction, imaging was performed weekly and imaging was performed daily during the last days of growth and first days of regression phases. During the late phases of tumor regression, we again reduced the imaging frequency to weekly or bi-weekly intervals.
Histology and H&E processing. The H&E staining refers to the mixture of Hematoxylin (Vector, H-3404) and Eosin Y (Fisher, LAMB/100-D), for nuclear and cytoplasm staining, respectively. Tissue histology and H&E staining were performed as outlined below, and following Cardiff et al. 65 . Tumors were extracted and fixed using 4% formalin (Sigma, HT501128), for 24 h at room temperature, followed by exposure to increasing concentrations of ethanol up until 100%, after which the sample is moved to xylene for clearing and finally embedded in paraffin. The paraffin block was then sectioned using a microtome to obtain single, 4.5 µm thick slices. Slides containing section slices were cleaned from paraffin using xylene (Roth, CN80.2), and then hydrated using a decreasing concentration of ethanol. The samples were exposed to Hematoxylin for 90 s and to Eosin for 30 s. After going through washes, the slides were once more put in contact with an ascending concentration of ethanol, ending on xylene. Slides were mounted using DPX mounting media (VWR, 360294H). As a consequence of the heterogeneous nature of the tumor phenotype, this staining protocol yielded stainings of variable outcome. In particular, only a small percentage of tumors had well-stained nuclei, which are a prerequisite to separate them from the cytoplasm and thus ensure accurate counting with commercial image processing software (HistoQuest software, Tissuegnostics, Vienna, Austria). This motivated us to devise a custom image processing pipeline (see Fig. 3a, b).
Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
The datasets generated and analyzed in the current study are in the supplementary data file or available from the corresponding author on reasonable request.