Robust Quantification of Regional Patterns of Migration in Three-Dimensional Cell Culture Models

Wound healing assays is a common two-dimensional migration model, with the spheroid assay three-dimensional migration model recently emerging as being more representative of in vivo migration behaviours. These models provide insight into the overall migration of cells in response to various factors such as biological, chemotactic and molecular agents. However, currently available analysis techniques for these assays fall short on providing quantifiable means to measure regional migration patterns, which is essential to allow a more robust assessment of drug treatments on cell migration in a chemotactic fashion. Therefore, this study aims to develop a finite element (FE) based pipeline that can objectively quantify regional migration patterns of cells. We have developed a novel FE based approach that is able to accurately measure changes in overall migration areas of 3D Glioblastoma Multiforme (GBM) spheroids that we generated using the primary cell lines from patients undergoing tumour resection surgery. We live-imaged the migration patterns of GBM spheroids and analysed them, first with the standard ImageJ method. We then performed the same analysis with the proposed FE method. When compared to the standard ImageJ method, our proposed method was able to measure the changes in a more quantitative and accurate manner. Furthermore, our regional migration analysis provided means to analyse the migration pattern seen in the phantom data and our experimental results. Our FE based method will be a a robust tool for analysing cell migration patterns of GBM and other migrating cells in various diseases and degenerations.


Introduction
Cell migration is an important process regulating physiological (e.g., capillary formation) and pathological (e.g., tumour spread) processes. A majority of in vitro migration studies use two-dimensional (2D) monolayer cell cultures assays.
Wound healing (scraping) and transwell migration (Boyden chamber) assays are some of the more common methods used to study in vitro cell migration [1,2]. More recently, the spheroid assay has emerged in the field of migration studies. This assay combines both two-2D and three-dimensional (3D) technologies to better represent the tissue morphology and cell-cell contacts without introducing cell-free regions [1].
Recent advancements of 3D cultures and experiments have allowed researchers to better mimic in vivo conditions such as tissue microstructures and cellular status, both of which can impact proliferation and migration [1,3]. The formation of spheroids is typically used to investigate tumour growth and viability and test anticancer agents [4][5][6]. Additionally, recent studies show that the spheroids structure provides a better representation of the physiological migration of the cells from a tissue-like structure, which is particularly beneficial in studying tumour cell migration, as solid tumours generally grow as a solid mass [7]. The spheroid assay begins with the attachment of the spheres on the extracellular matrix (ECM)-coated cell culture plate surfaces. Following attachment, cellular migration from the sphere can be influenced by various factors, such as serum [8,9], chemotactic agents [10,11], and migration inhibitors [12,13]. As a result, the cell movement can be measured microscopically and measured at different time-points. This assay can be applied to cell types that can form spheroids, such as neural stem cells [14,15], certain embryonic stem cells [16,17], induced pluripotent stem cells (iPSC) and their organoids [18][19][20][21], and certain solid tumour cells [7,22]. The migration of cells is quantified by measuring the area covered by the cells leaving the sphere [23] or by measuring the normalised radial migratory distance of the cells from the sphere relative to the original diameter of the sphere [24].
Cells in the wound healing assays are either stained or unstained and imaged at multiple time points for migration analyses. Analysis of wound healing migration assesses the overall migration of the cells to cover the cell-free region. In these cases, researchers typically use ImageJ or Metamorph analysis (Molecular Devices Corporation, Downingtown, PA) software to measure the distance between one side of the scrap and the other or area of the cell-free region [2,[25][26][27][28]. The spheroid assay is analysed in a similar way to the wound healing assay, as the spheres and cells are typically unstained and imaged at multiple time points for migration analyses. The collective migration of cells in this assay is quantified by measuring and calculating the area of cells that have migrated from the sphere [23]. Another way to quantify overall migration in this assay is to calculate the radial distance that the cells have migrated from the sphere and normalise it to the sphere's mean initial diameter [24].
Although the above methods are effective at assessing the overall migration of the cells radiating from a sphere, they fall short on providing a quantifiable regional migration analysis in which there are elements of directional migration. Regional migration analysis could prove advantageous as it provides greater insights into migration patterns, more robustly assesses drug treatments on cell migration in a chemotactic fashion, and provides more avenues of experimentation in 2D and 3D models.
One way of analysing regional migration is to create a computational model of the collective cell migration shape and measure the shape changes over time. In computational image analysis, shape changes can be quantified by measuring the amount of deformation between original undeformed and deformed states [29], which is often described as strain. This has advantages over other methods that measure quantitative shape features, such as area and length because the strain-based method allows more detailed and region-specific analysis in the cell migration patterns.
One method that allows such analysis is Finite Element (FE) models, which have been used extensively in various bioengineering related problems such as joint biomechanics [30,31], tissue deformation [32] and injury prediction [29]. FE analysis has also been used on a micro-scale, such as Saeed and Weihs, who simulated a 3D FE cell model and investigated the cell's morphology changes as it underwent uniformly applied compressions [33]. In our previous study, we have used FE analysis to investigate the deformation of cells [34] and cell-seeded 3D hydrogels and characterised the inhomogeneities in strain distribution [35]. Some also used the FE method in analysing cell migration, particularly during wound healing. Most notably, Zhao and colleagues developed an FE cell migration model that described individual and collective cell movement, as well as the change in cell morphology. The model was consequently used to examine the directionality and persistence of cell migration during re-epithelialisation [36]. However, there are no FE models developed to describe regional cellular migration from spheroid structures, including tumorspheres. Furthermore, many of these studies deployed FE analysis tailored to their specific purposes, which limits the applicability of those methods to other problems.
In this study, we live-imaged the migration 3D Glioblastoma Multiforme (GBM) tumorsphere cells and used an FE-based cell migration shape model to analyse various modes of migration readings. This method provides a robust method of regional migrational analysis, allowing researchers the ability to quantitatively assess migration patterns in non-uniform 2D and 3D migration models.

GBM Spheroid Formation
TICs were harvested for passaging and experimentation by first dissociating the cultures into single-cell suspensions using Accutase® (Invitrogen). In order to generate tumourspheres, TICs were plated at 5000 cells/well in a U-shaped sphere-forming 96-well plate (Thermo Fisher Scientific). Once plated, the cells were incubated at 37 °C with 95% air, 5% CO2 and incubated for 72 h to ensure proper cell adhesion and sphere formation prior to any experiments.
After incubation, TIC spheres required extraction from the U-shaped plates-The spheres were extracted from the well along with 50 µL of the media in the well into either a Matrigel-coated 12-well plate (for live imaging; 1-2 spheres per well). The plated spheres were then incubated at 37 °C with 95% air, 5% CO 2 for approximately 2-3 h to ensure successful sphere adherence to the Matrigel®-coated plate. Following adherence, new stem cell media was added into each well-the total volume per well in a 12-well plate should be 1 mL.

Live-Imaging of Tumoursphere Cell Migration
After the tumourspheres were plated, they were transferred to the temperature, and CO 2 -controlled Live Cell Microscope Incubator (Zeiss XL-3 stage incubation chamber) fitted with the Axiovert 200 M motorised inverted microscope and AxioCam MRm digital camera (Zeiss). The temperature and CO 2 levels were maintained at 37 °C and 5% respectively within the chamber and adequately humidified to prevent excessive evaporation. The microscope magnification was maintained at 10 × with 0.25 bright-field images using A-plan every 20 min using the AxioVision program associated with the imaging suite (Zeiss) over the time course of 48 h. Colour and exposure settings were set at Black and White settings to ensure a clear distinction and focus of migrating cells.

Migration Analysis of Tumoursphere Cells
After the images were captured, they were first imported into ImageJ. ImageJ migration area analysis began by highlighting the 'Find Edges' operation and thresholding. The area was subsequently measured using the magic wand tool and finding the net difference between the migration area and sphere area.
Additionally, ImageJ was used to draw boundaries around the migration area which was used to deform the mesh according to the boundary shapes. First, a reference mesh was generated using high order cubic Hermite basis functions, which preserve both the continuity of nodal values and their first derivatives. The resulting high order elements can describe smooth surfaces often found in biological structures with a lot fewer elements as each node has total 24 degrees of freedom [32,37]. Therefore the reference mesh, which had 32 elements and 72 nodes, has 1728 degrees of freedom. This mesh was then deformed to fit the extracted boundaries. This was achieved by orthogonally projecting a cloud of data points to the mesh surfaces, to which the reference mesh was deformed into using the least square fit algorithm (Eq. 1; Fig. 1) [37]. In this study, the drawn boundaries were 2D, so the z-dimension of the 3D mesh was fixed to ensure an accurate representation of the 2D boundary.
where z d are the coordinates of data points and 1d , 2d are the results of the orthogonal projection of data points to the surface of the mesh. F s is the penalty function which ensures the geometric smoothness of the mesh. The expression of the penalty function is given in Eq. 2: where i (i = 1…5) are penalty parameters. Each of these terms has a distinct effect on the final shape of the fitted mesh; 1 and 2 terms control arc-length, while 3 and 4 control the arc-curvature in the 1 and 2 directions, respectively. The last term ( 5 ) represents the face area. The resultant meshes were generated at different time points and imported for FE analysis. From this, the changes in the migrating areas can be obtained by solving the governing equation that describes the deformation of an object at different time points. This was achieved by solving the static Cauchy equation (Eq. 3).
where, ij are the Cauchy stress tensor components, x j are spatial coordinates. b j are body force components, and is the soft tissue density; both variables were neglected in this case.
To quantify the migration rates, displacement boundary conditions were used as the migration pattern of the cells was already known from the images at different time points. This prescribed the nodal positions at each time point and our finite element solver computed for the remaining degrees of freedom. Then the change in the migration pattern was analysed using strain, which was calculated in the form of Lagrangian strain using the theory of finite deformation elasticity as in our previous study [34,38]. The strain between original and deformed (i.e. time lapsed) meshes was computed using the deformation gradient tensor between the two, which is quantified by measuring the change in length of material segments [34,38]. Given a material vector in the undeformed state (dX), which is mapped to the deformed state (dx), the gradient deformation tensor (F) is defined as the following (Eq. 4) Strain in a deformed mesh is determined by measuring segment length changes, which is done by computing the square length (ds 2 ) for the deformed segment dx giving (Eq. 5) From this, we calculate the Cauchy-Green deformation tensor (Eq. 6): Finally, using the deformation, we derive the Lagrangian finite strain tensor (Eq. 7): The cell deformation and internal strain patterns in tumoursphere at each time point were described with the final Lagrangian strain tensor. As for the material properties, since our main focus is the analysis of the overall migration area, we ignored any nonlinear viscoelastic behaviour that might occur in the internal region of the cell aggregates. Therefore, the cell aggregates were modelled with a hyperelastic Mooney Rivlin material with a homogenous property. The parameter C 1 and C 2 values were obtained from the experimental study of the tumour tissue [39]. We used the maximum principal strain from the final deformed mesh to describe the migration rate. We used our in-house FE software called CMISS, which is freely available for academic use for migration analysis (www. openc miss. org).

Statistical Analysis
Graphing and statistical analysis were performed using GraphPad Prism 8 (GraphPad Software Inc.). Correlation analysis was achieved in GraphPad Prism 8. Graphs with error bars were displayed as means ± standard error of measurement (SEM), and statistical significance was set at p < 0.05. Two-way ANOVA with Dunnet's post-test group comparison was completed.

ImageJ vs FE-Based Area Analysis
The rate and the overall area of migration are key metrics analysed in cellular migration experiments. When involving the migration from spheroid structures, such as tumourspheres, ImageJ area calculation has often been employed to measure these parameters. The aim of this section was to validate the FE analysis capability to achieve the level of accuracy offered by these standard methods of analysis. In order to validate our analysis, the tumoursphere migration rate and overall area were assessed using both the established ImageJ protocols and our FE analysis protocol. Figure 2 shows ImageJ and FE area analysis of two migrating tumoursphere repeats. The correlation analysis between ImageJ and FE areas showed a very high Pearson correlation, with an R 2 of 0.9917 in tumoursphere repeat 1, and R 2 of 0.9835 in tumoursphere repeat 2 ( Fig. 2A-2 and B-2 respectively). This confirms the accuracy of area analysis done by our FE based method.

FE-Regional Migration Rate Analysis
Phantom data with an evident preferential migration was firstly used to show the capabilities of our FE based method to quantify the directionality of migration. The analysis method was further applied to a real tumoursphere sample with a more nuanced migration to evaluate its effectiveness in a real-life example. Figure 3 illustrates the preferential rightward migration of the phantom data, which is visualised by the red data points in Fig. 3A2. The data points were segmented into 60° quantifiable regions, with significantly greater migration distance detected in the 120° segment throughout the four time points (p < 0.0001) when compared to all other segments (Fig. 3B). The ability to delineate and measure regional migration rates are illustrated in Fig. 3C. For example, the rate at which the 120° segment migrated was significantly higher during the first timepoint (p < 0.0001), which faithfully represented the phantom data is analysed (Fig. 3A2). Overall, the phantom data results allowed for optimisation and validation of the FEbased approach to quantify regional migration.
In Fig. 4, the tumoursphere exhibited an unexpected protrusion of cells on the right side, which did not migrate as much as the rest of the tumoursphere throughout the first 24 h. Quantification of the tumoursphere migration showed a relatively low migration rate at the 120° segment (Fig. 4B), which corresponded to the lack of migration in the area of cellular protrusion. Moreover, the migration distance of the 120° segment was significantly lower than all other segments throughout the 24 h (Fig. 4B), while the migration rate was significantly lower than all other segments in the first 6 h (Fig. 4C). However, only the 120° segment showed a steady increase in migration rate over the first 12 h while all other segments decreased over the same period (Fig. 4C). Overall, our FE-based approach was able to detect and quantify the variations in a sphere shape and migrational patterns of cells from live tumourspheres.

Discussion
In this study, we have developed an FE-based method that can measure the area changes, as well as the directionalities of the change in a tumoursphere model. To our knowledge, this is the first study that used FEA in analysing 3D tumoursphere migrations. Our method was highly comparable to the area changes measured by ImageJ-a standard tool used by many investigators to assess the overall cellular migration. But more importantly, the use of the FEA method provides a number of advantages over the standard ImageJ based Fig. 2 Processing of a representative tumoursphere sample for ImageJ and FE model images. Area measurement using ImageJ was achieved by identifying the boundaries using the magic wand tool on thresh-olded images. FE area analysis required a Free Hand boundary region before it was imported into our FE analysis software as a mesh and the area was measured. White scale bar 250 µm method. First, the use of FE meshes allows the analysis of regional migration in GBM tumourspheres, as shown in our results with phantom datasets (Fig. 4). This type of regional analysis is not possible using ImageJ. Moreover, the use of FE strain as a parameter for migration rate provides a unique ability to quantitatively characterise directional changes in the migration patterns, which is not possible with simple area calculation as well (see Fig. 5).
Our pipeline is also able to section the tumoursphere into any desired segments during our analyses. This allows the investigation of migration from different radial regions of the sphere that could detect and quantify non-uniform migratory patterns that could have been missed using conventional overall area change models. Furthermore, the current pipeline involves the segmentation of raw images into boundary lines and the fitting of a mesh onto the boundary lines, therefore broadening the applications for this pipeline. For example, this pipeline can be applied to other migration models such as wound-healing/scratch assays [28,[40][41][42][43], as the boundary of migration can be defined, and a mesh can be constructed and overlayed. For instance, Zepecki observed that Lck-inhibition in tumoursphere cells resulted in reduced cellular migration when compared to the controls [44]. Additionally, Zepecki and colleagues analysed the tumour growth in a 3D brain reconstruction of a growing tumour within a mouse's brain using Neurolucida explorer (MBF Biosciences) [44]. Hypothetically, these observations could be analysed using our pipeline to provide a quantifiable measure of migration, both in vitro and in vivo. Moreover, all tumoursphere migration analyses from a recent literature survey revealed that only the overall migration area was used to quantify tumoursphere migration propensity [7,23,24]. Our pipeline can include additional analysis modes, such as overall migration rate and, more importantly, regional migration rate. As a result, this can provide a deep and potentially, novel insight into migration patterns of cells, especially in response to migration modifying agents, such as extracellular matrix and cell surface modifications [45][46][47][48], chemotactic agents [49][50][51][52] and small molecular agents [53][54][55][56]. Finally, the use of 3D FE model means that we can also 3D image datasets such as confocal microscopy images or magnetic resonance imaging in our analysis to characterise migration patterns in 3D, Fig. 3 Changes in area measured in ImageJ and FE analysis correlated very highly with each other. Overall area was measured using ImageJ and our FE based method over 24 h in 6 h increments on two migrating tumoursphere repeats (A1 and B1). Relative area for both ImageJ and CMISS applications were calculated by dividing the current area of the boundary by the initial boundary area at time = 0 h. The areas were correlated using the William Pearson Correlation analysis in Graphpad Prism 8 and the R 2 was measured (A2 and B2) which will provide a unique advantage that cannot be done with simple area calculations.
Although the study shows promise in analysing regional migration in spheroid models, there are areas where improvements are required. First, the assumption of our analysis is that the boundary motion is the major factor that influences the movement of the entire cellular aggregate. Hence it will not be able to describe the movements of the internal cells within the aggregates. This will require non-linear viscoelastic material descriptions. Yet parameters for such complex material descriptions for tumour cells are yet to be estimated, making it not possible to incorporate in the modelling framework. Moreover, our main aim is to characterise the overall migration pattern of the cell aggregates and our experimental data consistently showed that it is dominated by the movement of boundary cells. Hence our assumption was appropriate for the study purpose. Secondly, the current pipeline requires manual processing of the images, which hinders the efficiency and scalability of the method. Moreover, FE-based methods are not widely used to analyse cellular/molecular biological processes; hence it may be difficult for users without prior experience in FE analysis to apply our method to their studies. Finally, the multiple steps required to execute the entire pipeline has not been integrated into a single Fig. 4 Our FE-based approach accurately quantified regional migration of the phantom data. The phantom data mesh (A1) was utilised in our FE models to quantify the migration rate over time which was visualised in the migration map (Red indicates high migration rate while blue indicates minimal migration rate; A2). Data points in the migration map were segmented in 60˚ increments and averaged, before they were plotted and compared against each other (B). Every segment at each timepoint were compared against the 120˚ segment and statistical significance was assessed via Two-way ANOVA test with Dunnett's post-test; ****p < 0.0001 all segments versus 120˚ segment (B). Migration rate was calculated by finding the difference in migration distances between the current timepoint and the previous timepoint, and the difference was divided by the time that has passed between the two timepoints (C) streamlined process, hence making it more complicated for users without experience in scripting languages. Some of these issues are currently being addressed by uniting many of the pipeline's steps under one programming language such as Python since OpenCV and FE modules in Python are capable of accomplishing many if not all of the steps in the current pipeline. As a result, a program made under one programming language will be more accessible to the broader audience and will be more integrated than the existing pipeline. We aim to report this progress in our next publication.
In conclusion, we have developed a novel customisable FE-based analysis algorithm that can be to detect regional changes in cellular migration and accurately quantify various migratory matrices. This can pave the way for more automated analyses of brain tumour migration patterns both in in vitro experiments and in vivo medical images.

Declarations
Ethical Approval All experiments were conducted with ethical approval by the University of Auckland Human Participants Ethics and the Northern Regional Ethics Committee. The primary cell line used in the experiments were obtained with consent from patients undergoing tumour resection surgery at Auckland City Hospital. The images of the the migrating tumoursphere were taken every 6 h for 24 h (A1) and were analysed by our FE models to produce the migration maps (Red indicates high migration rate while blue indicates minimal migration rate; A2). The data points from the migration maps were plotted and graphed in segments to assess the regional migration of the tumoursphere. Every segment at each timepoint were compared against the 120˚ segment and sta-tistical significance was assessed via Two-way ANOVA test with Dunnett's post-test; **p < 0.005, ***p < 0.0005, ****p < 0.0001 all segments versus 120˚ segment, and different *p from different segment comparisons are outlined in the graph (B). Migration rate was calculated by finding the difference in migration distances between the current timepoint and the previous timepoint, and the difference was divided by the time that has passed between the two timepoints (C). White scale bar represents 250 µm Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.