4D imaging and analysis of multicellular tumour spheroid cell migration and invasion

Studying and characterising tumour cell migration is critical for a better understanding of disease progression and for assessing drug efficacy. Whilst tumour cell migration occurs fundamentally in 3 spatial dimensions (3D), because of the practicalities of growing 3D cell cultures and the complexities in measuring and tracking movement in 3D over time, most migration studies to date have performed analysis in 2D. Here we imaged live multicellular tumour spheroids with lightsheet fluorescence microscopy to determine cellular migration and invasion in 3D over time (4D). We focused on glioblastoma, which are aggressive brain tumours, where cell invasion into the surrounding normal brain remains a major challenge for clinical management. We developed a workflow for analysing complex 3D cell movement, taking into account migration within the spheroid as well as invasion into the surrounding matrix. This provided metrics that summarize the motion of the cells in different experimental conditions. As a proof of principle, we show how 3D tracking and analysis can be used to evaluate the efficacy of chemotherapeutics on invasion. The rich datasets open avenues for further studies on drug efficacy, microenvironment composition and stiffness, as well as collective cell migration and metastatic potential.


Introduction
The multicellular tumour spheroid is a well-characterised model in which cells are grown as three-dimensional (3D) clusters in the absence of any exogenous support. The spheroid model replicates important aspects of the tumour microenvironment such as gradients in oxygen, nutrients, pH, and cell proliferation and survival (1)(2)(3). The cellular and molecular heterogeneity of spheroids makes them a better model for solid tumours than homogeneous monolayer cultures. These pathophysiologically-representative 3D models are increasingly important in cancer research, especially for drug testing. The invasion of tumour cells into neighbouring tissue or distant organs is a hallmark of cancer and a major cause of mortality (4), which urgently needs to be tackled therapeutically. Good in vitro models for testing drugs that inhibit cell invasion are therefore essential. The majority of in vitro invasion assays to date, and associated drug screening, have been conducted using cells cultured as twodimensional (2D) monolayers on plastic and glass surfaces (5). These conditions fail to replicate important aspects of the cellular environment that occur in vivo, such as complex cellcell and cell-matrix interactions, matrix stiffness, and gradients in soluble factors such as oxygen, nutrients, growth factors, and cytokines (6). This explains why the results of in vitro drug screening often only poorly correlate with results obtained in vivo. Hence, drug testing in spheroids is emerging as a beneficial intermediate between testing in monolayer cultures and testing in animals (7). In previous attempts at measuring cellular invasion in 3D models, the spheroids were plated on, or embedded in, extracellular matrix (ECM) in 96-well plates and imaged using an inverted microscope (8)(9)(10)(11) providing information about the movement of some of the invading cells (in the x-y plane) but failing to capture the movements of cells within the spheroids. Imaging thick, living samples over long time courses can be challenging due to the effects of light diffraction, photobleaching and phototoxicity (12). Lightsheet fluorescence microscopy (LSFM), also known as single-plane illumination microscopy (SPIM), has overcome some of these challenges. By illuminating only the focal plane of the detection objective, the effects of photobleaching and phototoxicity are much reduced compared to confocal microscopy (13,14). LSFM also allows for very fast image acquisition from multiple angles, making it the ideal technique for imaging cell movement in a large, tightly packed spheroid, where z-stacks need to be acquired in rapid succession to facilitate the tracking of cells in space. Concurrently, the development of new software tools for spatial-temporal quantification of tumour spheroid dynamics (15) facilitates the analysis of such data rich 3D time lapse experiments. Here we used live LFSM of multicellular tumour spheroids to characterise cellular migration and invasion in 4D. We developed a workflow for analysing complex 3D cell trajectories with either commercial or open source software platform for tracking, taking into account migration within the spheroid as well as invasion into the surrounding extracellular matrix. This generated metrics characteristic of cell motion. These methods are applicable to a range of experimental conditions, including anti-cancer drug treatment. As proof of principle, spheroids were treated with 17-allylamino-17-demthoxygeldanamycin (17AAG), a heatshock protein 90 (HSP90) inhibitor, and cilengitide, an αvβ3 and αvβ5 integrin inhibitor with anti-tumour activity (16). Live LSFM combined with cell motion analysis allowed us to discern the effects of these molecules both on the cells invading away from the spheroid as well as on movement within the spheroid. Both effects may need to be considered P R E P R I N T to evaluate the mechanisms and efficacy of anticancer drugs in inhibiting tumour cell infiltration.

Methods
H2B-mRFP-U87MG stable cell line production. A nuclear marker is necessary to facilitate object recognition by the tracking software and enable automated cell tracking. We previously found that nuclear stains with dyes show poor penetration in spheroids, therefore, we created a cell line stably expressing a fluorescent histone H2B nuclear reporter protein by transducing U87-MG cells with pHIV-H2BmRFP. Lentiviral particles were generated by transfecting HEK293T cells with the following plasmids from Addgene (http://www.addgene.org): pHIV-H2BmRFP (Plasmid #18982), psPAX2 (Plasmid #12259) and pMD2.G (Plasmid #12260) in a 4:2:1 ratio. On day 1, HEK 293T cells were seeded in six 10 cm 2 dishes at a density of 1.5 x 10 6 per dish in 10 mL growth media (DMEM + 10% FBS + nonessential amino acids). On day 2, cells were transfected with the aforementioned plasmids using polyethylenimine (PEI) at a ratio of 2 µL: 1 µg PEI:DNA. On day 3, the media was changed and on day 5 supernatant was harvested, centrifuged at 500 x g for 10 min and filtered with a 0.45 µm PES filter before being subjected to ultracentrifugation following the protocol described in (17). 2.5 x 10 5 U87-MG cells were seeded in to a 25 cm 2 tissue culture flask and 250 µL virus was applied. Media was changed after 24 h. Spheroid culture and preparation. We have previously described a workflow for the growth, handling and imaging of spheroids made with fluorescently-labelled cells (18). The hanging drop method was used to grow compact spheroids from U87-H2BmRFP glioblastoma cells over 3 days (detailed protocol is available in (18)). Three day old spheroids were used to limit the effect of light scattering seen when imaging deeper (past~150 µm). At the beginning of imaging, spheroids were approximately 300 µm in diameter. Spheroids were prepared for imaging by aspirating the growth media and replacing it with 50% ice-cold Matrigel (Corning, Netherland) and 50% growth media plus 25 mM HEPES. For drug treated spheroids, 0.5 µM 17AAG (Abcam, UK) or 2 µM cilengitide (Adooq Bioscience, CA, USA) was added to the growth media. A Teflon plunger (Zeiss) was used to draw the spheroid into an optically-amenable fluorinated ethylene propylene (FEP) tube (S 1815-04, BOLA, Germany), which was sealed at one end using parafilm, mounted in the Lightsheet Z.1 (Zeiss) sample holder and then immersed in water at 37°C for the duration of the experiment.
Image acquisition. Samples were excited with a 561 nm laser and 10x excitation objectives; emitted light was collected through a 576-615 band pass filter using a 20x W Plan-Apochromat objective and 1x zoom. Images of the H2B-mRFP-labelled nuclei were acquired every 3 min with a zstep of 1.66 µm, starting approximately 50 µm in front of the spheroid. Lightsheet Z.1 Zen software (Zeiss) was used to acquire images using online dual side fusion with pivot scan. Images were detected using a pco.edge scientific complementary metal-oxide-semiconductor (sCMOS) camera.
Image pre-processing. Processing was done on a high-end workstation (cf. a high performance compute cluster), and thus reducing the size of the lightsheet data prior to analysis was essential due to the large file sizes produced (up to 0.75 TB for one experiment). A custom Fiji macro (see code availability below) allowed us to load one time point at a time and to down-sample the data by decreasing bit depth from 16 bit to 8 bit and by performing 2x binning (this is user configurable). A further processing step was carried out in Fiji prior to analyses in order to deal with sample drift (19), using the descriptor-based series registration plugin developed by Preibisch et al. (20). Registration was carried out using a translational transformation model and consecutive matching of images with a range for all-to-all matching of two.
Cell tracking. There are many excellent software packages available for tracking points over time (21). In this study, both the commercially-available Imaris 8.

Track ID to link features across frames
Although most software provides further analysis such as calculating overall track speed or intensity of spots over time, to give our analysis the widest utility, we wrote our processing scripts to assume no a priori track analysis information other than these five outputs. Imaris In Imaris, tracking was carried out using an autoregressive motion algorithm, with a maximum distance of 6 µm and a maximum gap size of zero. A filter was applied to remove any tracks that were shorter than eight time points (24 min). Once tracking was complete, the 'Position' spot statistics were exported into a CSV file. Trackmate The control datasets were also tracked using the Trackmate plugin (v 3.8.0) for Fiji, using the same parameters as above. The 'Spots in tracks statistics' output of Trackmate were exported into a CSV file. CSV outputs from Imaris or Trackmate can be imported into the MATLAB script for further processing.

Measurement of spheroid boundary and centroid.
The spatial location of the outer boundary of each spheroid was defined by taking the mean radial distance of the ten furthest out spots present in the first frame of the experiment and adding 5%. Because part of the analysis relies upon the isotropic nature of a spheroid, the coordinates of the spheroid centroid has to be provided by the user, which allows characterisation of the position and distribution of cells. This is requested from the user at the beginning of the data processing P R E P R I N T script. It is important to use the same software for estimating centroid position as for tracking, as different software may use different frames of reference to define the image origin. Using Imaris, three orthogonal clipping planes were created and positioned centrally to allow the read-out of the relevant position on each axis. The same result can be achieved in Fiji by measuring the X,Y & Z coordinates using the cursor position to read out on the status bar, or alternatively, using the centroid position from a segmented volume.
Analysis of tracking data. MATLAB (Mathworks, US) was used to process the tracking output. A processing script was used to calculate spot and track parameters and a second script was used to create a characteristic set of plots. Briefly, an output CSV file was selected and, when requested, the Cartesian coordinates of the spheroid centroid and the time calibration (sec / frame) were entered. The boundary of the spheroid was calculated as the mean radial distance of the outermost 10 features in the first frame of the movie (see section 8 of the processing script). As spheroid drift was corrected in image pre-processing, this was a reliable measure of the boundary (see Discussion section for commentary on suitability). One of the key functionalities of the processing script is to break apart individual tracks into sections that are within the original spheroid boundary (defined at the start of the experiment) and those outside the boundary. Track-level summary statistics are provided as the mean parameters (e.g. windowed spot speed, windowed straightness, etc.) for all inside and outside sections separately. The proportional time spent inside and outside can be inferred.
Plotting and visualisation. The main output of the processing script is a MATLAB workspace (.mat) file containing all of the calculated values for the features and tracks (see the code repository for a detailed description of the output tables). For meaningful visualisation, a second script was created, providing the summary plots shown in the figures. The plotting and output options of this script are controlled by a number of variables whose parameter space and functions are detailed in the script header (for example axis line widths, font sizes, etc). Other than the input and initialisation code section, each plot is written as a MATLAB section allowing the code to be run in completion or in sections to output only the desired figures. We have found Zegami (http://zegami.com) to be a useful tool for querying and investigating the large number of graphs that are produced by the script processing. The output supporting the datasets in this study are openly available and can be explored at http://zegami.liv.ac.uk under the project "Spheroid Tracking Data".
Code Availability. The pre-processing Fiji script and the processing and visualisation MATLAB scripts are available at https://bitbucket.org/davemason/ lsfm_scripts along with documentation and example tracking data. At time of publication the commit hash is 95a4329.

Results
Single track analysis in 3 dimensions. Previous attempts to describe the migration of spheroid-associated cells over time have largely focussed on gross morphological changes to the spheroid itself, using 2D imaging or 3D projections (23). Such depictions may miss important dynamics occurring at a single cell level or within the body of the spheroid. This is shown in Figure 1A, where a 2D (maximum intensity) projection of two time points is displayed. To characterise single-cell movement, we developed an analysis pipeline to quantify the migration of single cells through post-acquisition tracking and track analysis (see overlaid tracks in Figure 1A). In this experiment, we observed cell migration within the spheroid as well as invasion into the extracellular matrix (Supplemental Movie 1). We tracked cell movements and obtained tracks from 3 control spheroids. One such track, showing the movement of a single cell, is shown in Figure 1B. Yet, the challenge is to extract parameters that allow to capture the behaviour of the population of cells as a whole, in particular in the context of invasion. It is therefore critical to assess the distribution of single track or single spot data. We found both the speed ( Figure 1C) and straightness ( Figure 1D) to be informative statistics for describing the motion of cells within and outside the spheroid. Unlike track analysis in an embryo (24) or tissue, spheroids are isotropic. With a defined centroid, a boundary can be demarcated between the edge of the spheroid and the surrounding matrix (see Methods) and thus all tracks have a geometric frame of reference ( Figure 1E). Therefore, the track positions can be characterised by their radial distances, and sorted into two categories depending upon whether they remain within the spheroid boundary or cross the boundary and invade the surrounding matrix (see Figure 1F).

Cell Population Analysis .
One of the most basic features of metastatic cancer cells is invasion; i.e. their ability to leave their initial niche and migrate elsewhere within the close microenvironment. However, there is no clear metric to define and characterise invasion and the extent of invasion is often underestimated through the projection of a 3D+time experiment onto 2D+time either through the imaging (for example imaging a single z-plane) or post-acquisition using a projection in the z-axis. By locating cell nuclei, we captured an accurate picture of invasion into the surrounding matrix and represented this through a cumulative distribution plot (Figure 2A). This type of plot describes the location of the cells at the start of the experiment (blue line shows frame 1) and how the distribution of the cells radial positions changes over the course of the experiment. Compared to a 2D projection, one important advantage is that information in the z-plane is preserved. Two interesting insights are gleaned from this plot. Firstly, the increase in the plateau (Figure 2A Figure 2E, F). We can show, for example, a mean track speed (± standard deviation) of 27 ± 10 µm/h within the spheroid boundary (blue line Figure 2C) and 29 ± 12 µm/h in the surrounding matrix. Being able to query the data in this way is useful because many perturbations in migration, induced for example by a drug treatment or genetic manipulation, may have an effect as a function of distance (a possible proxy for cell density, cellcell contacts, effector molecule concentration gradient, etc.) or time (which may be relevant to the sensitivity of feedback loops or cell cycle dependence). In this model, distance is especially relevant given the composition and physical properties of the embedding gel. The normalised plots in Figure  2C, D are shown with two bins of spatial data corresponding to inside vs outside (based on the initial size of the spheroid -see details in Methods). They demonstrate that while there is no discernible difference in speed, the straightness is increased outside the spheroid boundary. The number of spatial bins are configurable by the user if, for example, when using a larger spheroid a greater number of bins would be more appropriate (see Methods section: Plotting and visualisation).

Analysis of radial components and track breakdown.
As spheroids are spatially isotropic systems, the movement of cells is best described in spherical coordinates in relation to the centre of the spheroid. Velocity vectors can then be calculated for azimuthal, elevational and radial components ( Figure 3A, B show examples of the different coordinate systems for a 2-dimensional system). The radial velocity was plotted as a function of the distance from the centre of the spheroid ( Figure 3C) or as a function of time ( Figure 3D). In the spheroid model, the sum of positive (i.e. outward) and negative (i.e. inward) velocities, tend to approximate zero. We wanted to investigate if the migration patterns were different for cells migrating within the spheroids or cells invading into the extra-spheroid space (or migrating back into the spheroid). The 'migratory' tracks (i.e. those which cross the boundary in either direction) into their inside and outside components were separated, taking into account multiple transitions. This allowed the comparison of all of the inside and outside components within an individual track. Despite the relatively few cells that migrate into the surrounding matrix over the course of the experiment (~15% of total number of cells detected), a slight positive skew can still be seen when the distribution of radial velocities is plotted separately for vectors inside and those outside the spheroid boundary ( Figure 3E). Furthermore, by pairing the components of single tracks, it is possible to get a clearer picture of whether individual cells behave differently in the matrix compared to within the spheroid. The distribution of the natural log of the paired ratios are plotted in Figure 3F. The data are presented this way such that zero (shown with a red vertical line) represents no difference when comparing inside versus outside speeds. The positively-skewed data (with mean of the absolute values being 6.6) indicate the absolute radial velocity of cells decreases when migrating through the matrix compared with their time inside the spheroid.

Testing anticancer drugs effects on cell migration and invasion.
To assess the suitability of this 3D single cell-tracking approach as a method to evaluate inhibitors of invasion, we treated spheroids with 17-allylamino-17demethoxygeldanamycin (17AAG) and cilengitide. Both drugs were added to the cell culture media during sample embedding in Matrigel. 17AAG has been shown to inhibit the growth and invasion of glioma cells by inhibiting heat shock protein 90 (HSP90), a chaperone protein that regulates numerous signalling proteins, including those involved in cell proliferation, survival, migration and invasion (25)(26)(27). Cilengitide is a cyclic peptide antagonist of integrins αVβ3 and αVβ5 that has recently been in clinical trials for recurrent and newly diagnosed glioblastoma (28). The evidence of the effect of cilengitide on invasion is mixed: while some reports document an anti-invasive effect (29,30), others have found that it acts to promote cell detachment, migration and invasion in 2D cell-based assays (31). Similar to Maurer et al, we observed that cells treated with cilengitide showed more rounding up, presumably due to the im-pairment of adhesion properties (Supplemental Figure 1A). To clarify the role of integrins in the different steps involved in invasion, we initially tested the effects of cilengitide on 2D cell migration. Cilengitide slightly increased the distribution of mean track speeds compared to untreated controls and triggered cell detachment from the cell culture dish (Supplemental Figure 1B). However, when 3-day spheroids were treated with the same concentration of cilengitide, invasion into the matrix was not inhibited ( Figure 4A, B and Supplemental Figure 2). The distribution of the windowed speed was comparable between control and cilengitide treated conditions ( Figure 4C and Supplemental Movie 2). The absence of cilengitide effects in the 3D model conflicts with the detachment and changes in migration observed in 2D cell cul-P R E P R I N T ture, and demonstrates the risks of extrapolating results observed in simple models to more complex ones, as the integrin composition is likely to change (32) (see discussion). Conversely, upon HSP90 inhibition with 17AAG, a lack of invasion was immediately evident from the time-lapse (compare Figure 4D to 4E after 16h and Supplemental Movie 3). The cumulative distribution graph ( Figure 4F and Supplemental Figure 2) was also markedly different than the control ( Figure 4A) and the cilengitide treated spheroids (Figure 4B) with only a small amount of curve shift at the very edge of the spheroid. Furthermore, the distributions of both the speed ( Figure 4G) and straightness ( Figure 4H) of the 17AAG-treated spheroids were decreased compared to control and cilengitide treated, suggesting both slower and more tortuous motion. In summary, we have developed a workflow for the analysis of spheroid cell migration and invasion data. The acquisition, post-processing and parameter extraction of single cell movement in spheroids allows effects of chemotherapeutic drugs on cell migration and invasion to be tested.

Discussion
While spheroids have previously been used to investigate invasion, these investigations have until recently been limited by the available imaging techniques and a lack of clear metrics by which to quantify and compare conditions. Projecting 3D imaging data or relying on 2D imaging alone can make the processing and analysis easier, yet this inherently decreases the accuracy of the measurements. We described here a workflow to aid in the acquisition and analysis of 4D data from lightsheet imaging, using spheroids derived from glioblastoma cells as a model system. By imaging and tracking individual cells, we were able to gain many of the broad insights provided in 2D acquisition systems as well as leveraging the detailed information gleaned from the 3D nature of the system. 3D cell microenvironment and migration speed. In line with in vivo observations, we recorded a high mobility of cells within gliobastoma spheroids. Under control conditions, the mean track speed of 27 µm/h within the spheroid boundary, correlates well with the results of Farin et al in which glioma cells injected into neonatal rat forebrains were tracked and evaluated to move at 25 µm/h in sectioned tissues when imaged (as with our experiments) every 3 minutes (33). However, the tissue microenvironment and its properties are likely to affect cell migration. Indeed, we found that track straightness was consistently increased in cells that were invading into the Matrigel compared with cells that were migrating within the spheroid ( Figure 2D). It has been previously described how the speed and mode of cellular invasion are modulated by properties of the ECM such as composition, matrix stiffness and pore size (34)(35)(36)(37). Our experimental and analytical workflow provides a precise and quantitative measurement of the cell migration properties that can be used to assess these changes. One question that remains unanswered is whether cells migrating out from a spheroid are creating new pores, or using existing holes created by cells that have previously left. With slight adaptations such as using fluorescently labelled matrix or cell body (see below), this workflow can be adapted to study this type of cell migration. Nevertheless, there are several advantages to nuclear tracking, chiefly that feature detection and segmentation during tracking is considerably more robust especially if cells are in contact. Touching cells expressing a soluble GFP or other cytosolic or membrane label are harder to segment within a complex 3D structure. Moreover, nuclear tracking allows the possibility of gathering additional information on cell division if required.

Drug effects on cell invasion versus cell migration.
To compare cell motion into the surrounding matrix to movement within the spheroid, a boundary has to be defined. The position of the boundary is somewhat arbitrary because there is a continuous transition in cell density between inside and outside of the spheroid. Here we opted for a simple definition of "inside" and "outside" based on the initial time point. We found this to be a robust and useful functional definition and furthermore can be used to reliably segment tracks as the "core" of the spheroid does not appreciably change in size during the experiment. Of note, the open availability of the code makes it simple for others to adapt the boundary calculation or indeed any other part of the data processing. Our results show that 17AAG prevented the invasion of cells into the Matrigel but interestingly did not prevent their migration within the spheroid, although the movement of cells inside 17AAG treated spheroids was slower and less directional than in control spheroids. 17AAG affects many signalling proteins through its inhibition of HSP90, making it difficult to pinpoint the mechanism through which invasion is inhibited (25,26,38). However, it has been shown that treatment with 17AAG and HSP90 inhibition modulates actin dynamics, thereby affecting cell migration and invasion (39). In contrast to 17AAG, we found that cilengitide was not effective in inhibiting invasion and, furthermore, we could not detect any differences between spheroids treated with cilengitide and control spheroids. Cilengitide reduces cell viability by causing cell detachment (anoikis). The sensitivity of U87 cells to cilengitide-induced detachment has been shown to be dependent on the composition of the ECM, with cells cultured on collagen-coated plates resisting detachment (31,40). Gritsenko and Friedl have recently demonstrated that β 1 and α v integrins represent the primary adhesion systems for glioma cell migration in different migration models (spheroids and organotypic brain slices) (32). However, when they compared the same integrin interference strategy in the different models, they observed a variation in the contribution of integrin-dependent glioma invasion. Their results suggest that there are alternative cell-matrix interactions in complex brain-like models, which drive different cell migration properties and require specific inhibitors. This illustrates the importance of replicating environmental factors, such as ECM composition, to obtain a valid assessment of P R E P R I N T anti-invasive or anti-tumour activity and points to the fact that even more complex models than single cell type spheroids might be necessary for in vitro drug testing. Nevertheless, the accurate measurement of cell velocity and movement, highly comparable to previous studies using glioma cells by Farin et al (33), supports the use of live spheroid imaging in drug screening, prior to costly investigations in preclinical animal models. Moreover, our approach has overcome some of the limitations faced by Farin et al, who were only able to track the movements of cells in the x-y plane (33). In conclusion, light sheet microscopy provides an excellent method for measuring cellular movements in multicellular tumour spheroids, providing not only high quality time-lapse images but also a framework to track and analyse these dense structures. Our findings with the two anti-cancer drugs tested indicate that it also holds promise for screening and investigating the mechanism of anti-invasive drugs, contributing to the principles of the replacement and reduction of the use of animals in experimental research.