1. Introduction

In recent years, in vitro spheroid and ex vivo organoid model systems have been developed to better mimic the physiological environments that exist in human diseases [1,2]. These three-dimensional structure (3DS) models allow for the accurate study of a range of in vivo physiological processes and tissue responses to mutation, damage and drugs. Such cellular models are considered significant technological innovations and are critical tools in cancer research and clinical applications such as anticancer drug discovery and personalized therapy [3,4].

Cell motion, either collective or individual, plays a pivotal role in many key biological processes such as embryogenesis, organ development, wound healing, cancer invasion and metastasis. Since it has been shown that changes in the spatial features of tumor architecture are fundamental markers of cancer growth and invasive behavior [5], the structural and dynamic information of 3DS tumor models is valuable.

In most studies that examine the invasive potential of tumors, one looks for dynamic cell motion of individual cells moving away from the 3DS into the surrounding matrix [6–8]. However, the internal motion of cells inside the 3DS is often overlooked. Although the movement of cells within spheroids has been shown to accompany the process of cells invading the surrounding matrix [9], the significance of this process as a precursor to tumor invasion has not yet been investigated. Accurate quantification of changes in cellular motility in 3DSs is critical in understanding cancer invasion. Thus, there is an urgent need for the detection, characterization and quantification of cellular motion taking place within a reacting 3DS.

The quantitative investigation of internal motion within 3DS has lately seen great advancements thanks to the development of a great variety of novel optical tools, among them white light interferometry, confocal heterodyne tomography [10] and phase-contrast microscopy [11–13]. These methods exploit the extraordinary sensitivity of small features in the coherent light scattering patterns to any movement of internally distributed scattering objects (i.e., cells) within the 3DS. Most, if not all of, the abovementioned methods and techniques require expensive sophisticated equipment which is nontrivial to operate, meaning they are, generally speaking, not practical for routine, wide-scale use.

In recent years, it has been suggested that the cellular mean velocity of internal cellular motion in 3DS objects may be extracted by classic dynamic light scattering (DLS) methods [14]. Moreover, an advanced multiple scattering MS-DLS scheme that exploits speckle dynamic analysis has been reported to yield quantitative velocity-distribution curves of cell populations within a biologically activated spheroid [15]. These results indicate that speckle analysis could be a practical, affordable investigational tool for the study of spheroid/organoid reaction.

However, the more generalized elaborate mathematical formulations such as those presented by Brunel et al. [15], may obscure some hidden relations due to their relative complexity when considering multiple scattering processes. Therefore, we suggest an analytical approach for the investigation of these same speckle patterns, with a significantly simplified approach towards data extraction and interpretation. According to our hypothesis, under the limitation of a small angle, single-scattering model, it would be suggested to directly extract various motion parameters, indicative of the uniformity of cellular migration directionality and the degree of order contained within apparently random cellular motion. Under appropriate experimental limitations we present a relatively simple algorithm that may enable a faster, on-line characterization of 3DS internal dynamics, with a potential impact in areas beyond the scope of this study, including biological and medical research.

1.1 Theoretical background:

1.1.1 Laser scattering pattern generated by a 3DS built of a large population of biological cells

Laser scattering speckle pattern is the general term for the optical phenomenon in which an irregular, granular pattern is generated by the interaction of a coherent illuminating beam with an irregular-shaped specimen. In the same manner, when laser light impinges on a large number N of randomly located biological cells within a 3DS (Fig. 1), the angular distribution of the scattered electric field is given by [16]:

 

Fig. 1. Formation of laser scattering speckle pattern. A coherent laser beam of wave-vector ${\bar{K}_{in}}$ impinges on a 3DS located in the origin of a Cartesian coordinate system. The scattered wave-vector ${\bar{K}_s}$ is defined by scattering angles ${{\boldsymbol \theta }_{\boldsymbol x}},{{\boldsymbol \theta }_{\boldsymbol y}}{\boldsymbol \; }$ at x and y axes correspondingly and total scattering angle ${{\boldsymbol \theta }_{\boldsymbol s}}$. For each cell within the 3DS at location ${\vec{r}_q}$ there exists a corresponding phase ${\overrightarrow {\delta k} _q} \cdot {\vec{r}_q}$.

Download Full Size | PDF

In the following, we focus attention on small scattering angles regimes where one may expect multiple scattering processes to have only minor effect on the speckle pattern. Accordingly, one may approximate $sin\theta \to \theta ,\; cos\theta \cong const.$, thus Eq. (3) becomes:

As mentioned above, the coordinates ${x_q}$ of a given cell, are generally random. Thus, the distribution of the electric field itself becomes a statistical entity. To describe the electric field auto-correlation function designated by $A{C_{EF\; }}({\delta {\sigma_s}} ),\; $ we define:

It can be shown that for a large number $N \gg 1\; $ of phasors then Eq. (5) may be approximated by the following [17]:

The phase difference $\delta \varphi $ in radians is assumed to be smaller than unity, which is to say that the expression obtained for the correlation function is valid in the vicinity of its peak.

The approximation presented in Eq. (6) may be exploited to generate some fundamental expressions that will be used in the next sections. The first, is the spatial (or angular) autocorrelation functions $A{C_{EF\; }}({\delta \theta } )\; $ of the scattering pattern created by a 3DS which is composed of a large number N of randomly located cells but fixed within it in time. The second to be introduced, $A{C_{EF\; }}(\tau )$, is the temporal autocorrelation function of the time varying speckle pattern generated by a 3DS composed of cells undergoing diffusive or non-directional motion (which is closely related to the– DLS experimental approach to be discussed later).

To calculate the angular autocorrelation $A{C_{EF\; }}({\delta \theta } )$, we first use Eq. (2) to find the average phase difference $\delta \varphi ({\delta \theta } )$:

The expression in Eq. (8) yields the well-known [17] correlation angle (defined as $\delta \theta $ that decreases $A{C_{EF\; }}$ to ${e^{ - 1}}$):

In random, non-directional motion, ${\tilde{x}_q}(t )$, the average displacement of a typical cell is given by [18]:

Again, with the well-known [18] temporal autocorrelation decay rate constant, ${\gamma _{rand}}$ is given as:

1.1.2 Representation of synchronized motion cells within a 3D biological cell object (3DS)

Invasion is a fundamental step in tumor progression toward metastasis. We assume that the motion of cells within a spheroid, precedes their movement into the surrounding matrix. Cellular motion inside the cancerous MCS might be expected to be synchronized (or coordinated), where cells start to invade into nearby ECM, upon a specific stimulus. A simple model of 3DS inflation or “stretching” in the x direction may be geometrically defined by the linear transformation:

Again, the auto-correlation function may be found by applying Eq. (8). The phase-difference this time, should be found using the full differential of ${\varphi _q}({\theta ,t} )$, namely:

It should be noted that the differential in Eq. (15) depends both on time t and scattering angle ${\theta _s}$. The physical meaning of the above, is that the temporal and angular correlation decay parameters are correspondingly time and angle dependent. We therefore define a ‘window’ or an experimental domain at which we collect the speckle data:

With the above assumptions, the relevant phase difference becomes

1.1.3 Composed internal motion: synchronized and random motion inside a 3DS

To further adapt our model for a more realistic biological 3DS, we incorporate the possibility of a random component occurring in the motion of cells comprising the 3DS. When cellular motion consists of both synchronized, inflation motion and random or diffusive components, we may write the time varying phase scattered by a cell q, as follows:

The expression obtained above is what we call the Spatial-Temporal Auto-Correlation (STAC) function for a 3DS composed of cells with both random and invasive components of motion. The importance of this result is in the fact that it enables a generalization of the conventional DLS approach, thus providing much more information regarding the 3DS under investigation.

To represent DLS experimental scenario in our novel STAC formalism, we do the following: Taking the case of non-synchronized, purely random cellular motion we apply Eq. (22) by setting ${\alpha _{x - \; axis}} = 0$ to obtain:

The corresponding STAC distribution map of (18) is shown in Fig. 2(A). Gaussian angular (horizontal) cross-section is seen with an exponential decay at the temporal (vertical) cross-section axis. Any other cross-section at lines parallel to the axes will produce the same Gaussian or exponential decay function, all centered at the origin.

 

Fig. 2. Spatial-Temporal correlation function distributions (STAC) with cross-sections at different points of angle and time. A. for the case of a random spatial and temporal motion of cellular population the STAC distribution shows Gaussian cross-sections in the angular direction and exponential at the temporal cross-sections all centered around the main axes. B. STAC distribution for the case of both random and synchronized motion of contraction along the measured axis.

Download Full Size | PDF

The image in 2A indicates that only at the spatial-temporal vicinity of a typical speckle (=origin), one would expect to find a significant trace of that same speckle. However, as we move further away from that typical speckle (on both space and time axes), it is expected to fade. As mentioned, this behavior corresponds to a standard DLS experimental scenario, where the main parameter required to describe the intensity fluctuations (presumably, related to vitality of the illuminated cells in the 3DS) are speckles typical size $\delta \theta _{speckle}^{}$ and their correlation (=fading) time $\tau _c^{}.$ However, the STAC function we get when setting a synchronized motion component in Eq. (22) so that ${\alpha _x} \ne 0$ stands in significant contrast. The two-dimensional distribution obtained, tends to spread at directions not parallel (to the time and space axes. As seen in Fig. 3(B), different parallel cross-sections yield different correlation curves which are not centered at the origin but are shifted by both temporal and spatial displacements.

 

Fig. 3. cross-sections of the STAC functions presented in Fig. 2(A) and 2B, illustrate DLS decay time as a special case of STAC formalism. A. left panel: with pure random motion of particles, a standard DLS correlation function decays, so that at different times (${t_0} - {t_3})$ one measures exponentially decaying smaller values of $A{C_{DLS\; }}(\tau )$. This situation is formally expressed as a fading of the peak of the STAC distribution, at the origin of angular correlation axes $\delta \theta = 0$. A. right panel: the corresponding sample values at $\delta \theta = {0_{}}$ at consecutive times are seen, which form the functional dependence of $A{C_{DLS\; }}$. B. left panel: in a spatially random but also temporally synchronized motion, then temporal cross-sections (${t_0} - {t_3})$ samples of the spatial-temporal function yield a drift of STAC peak. B. right panel: the samples at the $\delta \theta = 0$, again, represent the DLS curve with decay time ${\tau _{DLS}}$ that evidently shorter than the actual, ${\tau _{STAC}}\; of\; $ moving STAC correlation peak.

Download Full Size | PDF

Further discussing the implications of the above, Fig. 3(A), presents a plot which corresponds to the plot in Fig. 2(A). In Fig. 3(A), a purely random diffusion process takes place, with angular cross-sections shown for a set of evenly spaced correlation delay times (${\tau _0} - {\tau _3}$). We then get a set of Gaussian distributions, all centered (as implied in Fig. 1(A)) at the origin, but with decaying amplitudes. Decaying values of the corresponding DLS correlation function are obtained by sampling these Gaussian functions at $\delta \theta = 0$, for the set of incremental delay times (${\tau _0} - {\tau _3}$). By doing this, an exponential DLS function is obtained as seen in Fig. 3(A) right panel. On the other hand, in the presence of synchronized motion, as shown in Fig. 3(B) left panel (corresponding to Fig. 2(B)), the Gaussian cross-section distributions for the same set of delay times are gradually shifted so that their DLS value sampled at $\delta \theta = 0\; $ might be very different from the STAC correlation value expressed by their peak value of gradually shifted Gaussian functions.

As seen in Fig. 3, STAC analysis may be used to differ random motion from synchronized. In the right panels of 3A and 3B, are shown the cross-sections of the spatial-temporal correlation distribution of random and synchronized motions, correspondingly. As shown in the Fig. 3(B) right panel, DLS time-correlation function $\textrm{A}{\textrm{C}_{\textrm{DLS}}}(\mathrm{\tau } )$ is obtained by sampling at the x = 0 while $\textrm{A}{\textrm{C}_{\textrm{STAC}}}(\mathrm{\tau } )$, follows the peak as it moves. Therefore, whereas the apparent DLS auto-correlation time may look similar to that seen in the Fig. 3(A) right, neglecting the presence of order in the motion. $\textrm{A}{\textrm{C}_{\textrm{STAC}}}(\mathrm{\tau } )$ indicates the difference between the random and ordered motion.

It is thus demonstrated, that when inspecting only DLS correlation function, then relevant information could potentially be overlooked. We would like to show how such valuable biological information can be extracted by exploiting our novel STAC function approach.

1.1.4 Extraction of analytical expressions for DLS and STAC correlation decay times

According to the above explanation, we may now use Eq. (22) to find correlation time ${\tau _{DLS}}$ for the general case of random and synchronized cell motion as the time it takes the function $A{C_{\sigma \tau \; }}({\delta \theta = 0,\tau } )$ to decay to ${e^{ - 1}}$:

Spatial-temporal correlation time (to first approximation with $\gamma _{synch}^{} \le {\raise0.7ex\hbox{${\gamma _{rand}^{}}$} \!\mathord{/ {\vphantom {{\gamma_{rand}^{}} 2}} }\!\lower0.7ex\hbox{$2$}}$) is then given by:

To get the decay time ${\tau _{STAC}}$ of the suggested $\textrm{}A{C_{STAC\textrm{}}}(\tau )$, we note that with STAC experimental approach we trace the moving peak (or maximum) of the correlation function. Therefore, as illustrated above in Fig. 3(A) and Fig. 4, at each point of time the following prevails:

 

Fig. 4. Optical set-up and related images: A. The laser beam illuminates a specimen located on the object plane. Two cameras (Cam 1 and Cam 2) record the events in the object and its speckle pattern, correspondingly. B. physical model for 3DS is composed of about 100 polystyrene beads with 6 µm diameter, that are confined inside an oil drop. Laser beam of about 300micron diameter illuminates the whole structure of the 3DS phantom to produce a speckle pattern. C. the corresponding speckle pattern generated by laser illumination of the 3DS, to be analyzed. The Range of Interest (ROI) is defined by the constraints in Eq. (16).

Download Full Size | PDF

By doing so, the exponent at the right of Eq. (22) vanishes, leaving us with only:

Therefore, in such case, ${\tau _{STAC}}$ becomes:

At this point, it seems reasonable to define a parameter which characterizes the extent of synchronization of the motion within the specimen, this Motion Synchronization Parameter (MSP) may be defined by:

The extent of synchrony of motion, is therefore extracted by finding both ${\tau _{DLS}}\; $ and ${\tau _{STAC}}$. The other important parameter is the rate of invasion-inflation ${\alpha _{x - axis}}$. This parameter is found by Eq. (26) to be:

The meaning of Eq. (31) is that the fundamental signal of a STAC-based measurement is the size-change increments rather than size. However, these size increments may be integrated to yield a measure for the extent of overall inflation. These aspects will be discussed more thoroughly in the experimental result section.

2. Materials and methods

2.1 Spheroid culture procedures

Generation and culturing spheroids: 3D multicellular models of cancer spheroids will be used to validate the STAC theory. Spheroids derived from cervical cancer cell-lines were grown in the wells of the microtiter imaging plate using an unscaffolded 3DS cellular system as previously described [7,19–20]. Tumor cells were seeded in the wells of the agarose-based microtiter plate and allowed to spontaneously create spheroids. During the first day after seeding, individual cells aggregate and form multicellular clusters, which continue to grow over the following days, each in an individual well. For analysis of invasion capacity, mature tumor spheroids (48 h in culture) on the imaging plate were embedded in an extracellular matrix composed of type I rat tail collagen as described [7,20–21].

2.2 Experimental extraction of a 3DS and simulation model

Both numerical simulation tools and experimental setup are used to demonstrate practical implications of the STAC approach proposed above. All simulations presented hereafter were produced by means of a Beam Propagation Method (BPM) based code, which enabled modeling a variety of experimental arrangements and in-spheroid motion, hence producing some dynamic speckle patterns that could be used for assessment and optimization of possible data extraction algorithms. The experimental setup in Fig. 4 is designed to satisfy the scenario described in the theoretical section. Figure 4(A), presents the optical layout with a laser source (Melles Griot 05 LPL 343-1 Helium Neon Laser Supply with Hughes 3222H-PC Neon Laser) aligned to illuminate the specimen by applying a Gaussian laser spot with a waist just larger than the diameter of the object under investigation. This object could be either a calibration phantom or a living spheroid. The specimen shown in Fig. 4(B) is a phantom that enables mimicking a dynamically coordinated motion within a spheroid. This phantom was composed of about one hundred 6 µm diameter polystyrene beads (Molecular Probes, USA) arrested in an ∼120 µm oil droplet. By applying certain mechanical manipulations to the specimen, internal motion could be induced within the phantom. Images of the 3DS model could be recorded by means of a camera (Fig. 4(b): Cam1). Illuminating of the above-described 3DS model, creates a speckle pattern which was collected and recorded by the second camera (Fig. 4(c): Cam2). According to the considerations in Eq. (16), speckle pattern images were collected during a defined time period $\Delta {t_{collect}} \ll \alpha \textrm{ }_x^{ - 1}$ and within a defined angular zone ROI (in Fig. 4(C)) $\Delta {\theta _{collect}} \ll {\theta _s}$. As we shall next see, this procedure was applied on both actual experimental results and those produced by the numerically simulated model.

2.3 Numerical model:

The numerical model was designed both for the purpose of validating physical and optical assumptions made along the theoretical part of the work, and for optimizing data analysis procedures to be executed on the actual experimental data. The core of the numerical analysis was based on a 3DS model usually composed of several hundreds, to a couple of thousands of cells, any of which may have any size and shape, which is assumed to be (but not necessarily) constant, during motion. The model used an angular spectrum beam propagation method to propagate a Gaussian illuminating beam through the whole sample as illustrated in Fig. 5. The three-dimensional numerical model of the spheroid had a lateral resolution of 0.5micron and 2micron resolution on the beam propagation axis. A detailed numerical representation of the experimental set-up was done by taking into account the illumination beam size and wave front Fig. 5(A) and by representing the 3DS by its cell locations, sizes and shapes. Beam-propagation through a 3DS was applied by a split-step method, wherein layers were set 2 µm from each other. The typical diameter of spheroids in our numerical models was about 300 µm with thicknesses of 100-200 µm, cell size in the range 10-14 µm of which calculated transport mean free path ${l^\mathrm{\ast }}\textrm{}$ was at the range 30-40 µm. Cellular motion within these spheroids could be totally random, ballistic or partially random. By propagating the illumination beam through the described thick model, the effect of multiple scattering process could be measured at any angular location of the speckle pattern by comparing the speckle dynamics of an equivalent thin spheroid of the same structure, except for its thickness. Out-coming wave front (pale blue arrow) was then propagated through a lens and a beam-stop (Fig. 5(C)) and guided into the receiving plane at Fig. 5(D). Then, the cells inside the 3DS numerical model are induced to move according to some pre-selected motion-rule, such as motion related to viability (random motion) or to radial “ballistic” invasion (radially oriented inflation or axial inflation). Then, for any such scenario, speckle pattern sequence is recorded to be processed later by any suggested algorithm, be it “classic” DLS, STAC or another method.

 

Fig. 5. A. In-coming laser beam is numerically modeled and set to properly illuminate the area in which a B. 3DS spheroid model located on the specimen plane. The scattered field distribution is then collected by C. system optics composed of a lens and a beam-stop which guided the light into the camera sensitive receiving plane D. by which fully developed speckle pattern may be recorded and stored as stack of E. consecutive speckle-patterns to be analyzed by motion-extraction algorithm.

Download Full Size | PDF

In each of the 3DS configurations, a different speckle pattern is created. The numerical model enabled investigation of the extent of correlation signal to be expected in various cases such as spheroid cells size, shape homogeneity and thickness of 3DS, while as mentioned, multiple scattering effect is also taken into account. The numerical model was adapted to yield a speckle pattern sample that matches the field-of-view collected by the experimental setup optics.

3. Results

In order to assess the validity of Eq. (19) in the theoretical section, presented in Fig. 6 are results of numerical simulations for two cases: a purely random motion and a hybrid motion, namely, with both synchronic and random components. Spheroid model had a diameter of 300 µm and thickness of 200 µm. We have found through inspection, that for the 3DS simulated model described above, then at regime ${\boldsymbol \theta }_{\boldsymbol s}^{} \le \textrm{}{35^o}\; ,$ there was less than 5% error in the values of extracted parameters ${\alpha _x}$, ${\tau _{STAC}}$ and $\textrm{MSP}$, due to multiple scattering.

 

Fig. 6. A. spatial-temporal speckle pattern (STSP) for purely random motion. B. STAC distribution generated from A, that has only one meaningful feature, namely the decay along temporal axis ${\tau _{DLS}}$ C. STSP for partially synchronized and random motion. D. STAC distribution generated from C,with added features to those in B. These features indicate decay times along both the spatial and temporal axes. Added information contained in STAC decay time ${\tau _{STAC}}$ and the slant angle ${\alpha _x}\; $ may be exploited for a richer characterizing of internal motion in 3DS.

Download Full Size | PDF

It should be noted that the STAC functions obtained by the numerical simulations as presented in Fig. 6(c) and (d) closely resemble those predicted theoretical distributions presented in Fig. 2(B). However, the above results are still qualitative. For the demonstration of quantitative results we turn to our experimental system and the polystyrene bead model. The Spatial-Temporal speckle pattern (STSP) for the said two cases, obtained of speckle images contained in small scattering angles are presented in Figs. 6(A) and 6Cthe. They are generated by building a three-dimensional array of data having two angular coordinates and one temporal coordinate. The STSP such presented in Figs. 6(A) and 6(C) are only cross-sections of the mentioned arrays along a selected angular axis, done only for clarity of presentation. The corresponding auto-correlation functions (STAC) were produced and are presented in Figs. 6(B) and 6(D), correspondingly. Statistical and dynamic information may be extracted from these STAC distributions. Specifically, note the slant of the STAC function and the correlation times ${\tau _{DLS}}$ and ${\tau _{STAC}}$ in Fig. 6(D). According to Eq. (30) and Eq. (31), these may be further used to draw a quantitative Motion Synchronization Parameter (MSP) evaluation and the quantitative effective inflation rate ${\alpha _{axis}}$, within the investigated object.

As stated above, an identical procedure was applied for both actual experimental data (Figs. 7(A), 7(B) and 7C) and simulated model (Figs. 7(D)., 7E. and 7F.). Given a sequence of events in the investigated 2D object (Fig. 7(A)), a corresponding sequence of speckle images was stored in a 3D array (Figs. 7(B) and 7(C)) in which the cross-section along a preselected axis of consecutive speckle images forms the Spatial-Temporal speckle pattern (STSP) seen in Figs. 7(C) and 7(F).

 

Fig. 7. Top Row: Extraction of spatial-temporal speckle patterns from a set of images for the experimentally squeezed 3DS and its numerical model. A. image of the 3DS object under investigation as recorded by Cam1. B. sequence of speckle-pattern images as recorded by Cam2. C. Spatial-Temporal speckle pattern (SPST) generated from B, is built by accumulation of intensity cross-sections along $\theta $ axis (as shown in B) for different times. Bottom Row: Corresponding images obtained by numerical model. D. numerical model of 3DS object presented in Top Row. E. sequence of simulated speckle-pattern images. D. simulated Spatial-Temporal speckle pattern generated from simulation E.

Download Full Size | PDF

The final step in the procedure would be to generate the ST auto-correlation function. The standard method for applying such operation on an image ${I_{ST}}\; $ is via Fourier transform, namely:

The resulting distributions $A{C_{ST}}$ or STAC such as in Fig. 6(B) and 6(D) obtained by the above, are in fact our proposed novel Spatial-Temporal Auto-Correlation (STAC) distributions.

As described earlier, the objects presented in Fig. 7(A) were composed of about 100 polystyrene beads. These were arrested in an oil droplet with a diameter of about 120µm. Being kept under a microscope cover-glass with a constant volume, its thickness could be slightly perturbed by gentle pulses of pressure applied on the cover glass. In reaction to each such pulse, the beads would move, each with a given displacement (amplitude and direction) and then, after the pressure is released, return very close to their original location.

A very important aspect of the internal motion described above, was that the effective size of the bead cluster was changing with time. This may be observed in Fig.8b, where beads to the right of the phantom had a larger displacement than the beads to its left. The corresponding speckle pattern is expected to follow the same rational and present some typical directionality in the spatial-temporal correlations.

The whole process was recorded both by Cam1 and Cam2. Speckle dynamics were recorded by Cam2 for later processing stages, whereas Cam1, collected dynamic and geometric data regarding microscopic specimen, for generating a full numerical model of specimen at image plane. Then, by means of a beam propagation method, it was possible to reproduce a simulated speckle pattern based on the whole optical layout of the experimental system as depicted in Fig. 5. Figure 8, shows the relation between the physical model (Fig. 8(A)), a visualization of the numerical chart for the model (Fig. 8(B)), which shows locations of the different beads inside the model and the displacement vectors attached to each bead, and the emerging simulated speckle pattern with designated axis which corresponds to the directionality of beads in the object plane (Fig. 8(C)).

 

Fig. 8. 3DS physical model and visualization of its numerical motion chart. A. an example for a typical microscopic image of the 3DS internal arrangement of the beads composing it. By tracking the locations of these beads during the pressure pulses applied on the oil-drop then. B. motion chart is constructed of the polystyrene beads that enables implementation into the numerical model. Next, from this chart, a highly accurate simulation model was generated, that produced. C. simulated speckle patterns dynamics. These patterns can then be subjected to the STAC procedure as depicted in Fig. 7

Download Full Size | PDF

The spatial information as seen in Fig. 8(B), was then implemented into the numerical model to produce STAC speckle data, to which STAC function was calculated.

In Fig. 9, all stages of the data processing method are shown for both experimental results and synthetically produced data, obtained by the numerical model described above. Figure 9, panels a and d show the STAC speckle images for experimental and simulation correspondingly, to each of the latter, corresponding STAC function is calculated (Fig. 9, panels B. and E.) from which auto-correlation curves $A{C_{ST}}$ are deduced (Fig. 9, panels C and F). Specimen dynamic parameters Motion Synchronization Parameter (MSP) and effective inflation rate may be calculated from times ${\tau _{DLS}}$ and ${\tau _{stac}}$ and ${\alpha _{axis}}$.

 

Fig. 9. Experimental and simulated data analysis procedures: A. and D. are experimental and simulated spatial-temporal (ST) speckle images correspondingly, for which STAC distribution B. and E. were calculated. The corresponding $A{C_{ST}}\; $ experimental and simulated auto-correlation curves are presented in C. and F. As expected, due to synchronized motion, both experimental and simulated $A{C_{ST}}$ curves show deviation between DLS-type and STAC-type correlation decay times.

Download Full Size | PDF

The values of MSP and inflation rate ${\alpha _{axis}}\; $ obtained by the above procedures, may now be compared quantitatively. As seen in Fig. 9, panels C and F, the various $A{C_{ST}}\; $ curves have different decay rates, as predicted theoretically above. After plotting $A{C_{ST}}$ on a semi-logarithmic scale, ${\tau _{DLS}}$ and ${\tau _{stac}}$ and ${\alpha _{axis}}$ are extracted and all necessary values for a quantitative examination of our method are available.

3.1 Tracing incremental structural changes by tracking of correlation peaks

An alternative and somewhat more practical approach in processing speckle data of the type presented above, would be by tracing the cross-correlation peak trajectory or correlation peak “drift” that was denoted by $\delta {\theta _{peak}}\; $ in Eq. (31). The effective size increments of the 3DS under investigation may be directly obtained by continuously monitoring the time dependence $\delta {\theta _{peak}}(t )\; $ of the correlation peak dynamics. Thus, by applying cross-correlation operation on consecutive speckle images obtained along the “cellular” motion, one may define a motion increment related to structural changes taking place inside the sample. Thus, a two-dimensional “trace” may be produced for the trajectory of the speckle correlation peak made over time, indicating possible inflation-contraction evolution. Figure 10 presents a trace recorded for one cycle of pressure-relief of the phantom.

 

Fig. 10. Tracing of speckle correlation peak motion obtained by squeezing the phantom. A. a two-dimensional trace of speckle correlation peak dynamics shows that most motion occurred along the X-axis (which is the expected axis according to Fig. 8(B)). The quantity of motion is given in terms of accumulated differential peak displacement, integrated over the whole period of interest. B the time series correlation peak motion along the x axis.

Download Full Size | PDF

 

Fig. 11. Spheroids that were subjected to A. control, and B. hypotonic shock solutions.

Download Full Size | PDF

The information contained in such traces as depicted in Fig. 10, enables sensitive, fast and meaningful assessment of ‘cellular’ dynamics within the phantom. For example, in the case of the phantom presented in the previous section, the time scale of the internal motion was about a hundredth of a second. Nevertheless, the ongoing internal motion could rapidly be detected and interpreted with by our simple correlation-peak tracking algorithm. As shown in Fig. 8, the mean displacement of the beads in the phantom was about 2-4 µm. By applying Eq. (31) to the dynamic speckle pattern data, the results presented in Fig. 10(B) show a displacement trace reasonably close to findings obtained by visual inspection of image data. Instantaneous mean velocity may then also be calculated by the rate of change in peak location ${\raise0.7ex\hbox{${\delta {\theta _{peak}}(t )}$} \!\mathord{/ {\vphantom {{\delta {\theta_{peak}}(t )} {\delta t}}} }\!\lower0.7ex\hbox{${\delta t}$}}.$ The evident potential of this method is in its sensitivity to internal motion, eliminating the need for special image-analysis hardware, large optical layout.

3.2 Spheroid osmolarity

Cellular responses to changes in the osmolarity of the cell’s surroundings have been widely studied [20]. These effects have also been measured in 3DS objects [21].

In the final section of this paper, we introduce the results obtained for spheroids which were subjected to hypotonic osmotic pressure, as depicted in Fig. 11. In principle, such an operation should induce a swelling process of the cells in the spheroid, thereby causing the spheroid itself to swell. Spheroids were placed on a microscope plate and subjected to either hypotonic (double distilled water – DDW) or control solutions (isotonic 280 mOsmol). Then, spheroids were investigated by both microscopic image analysis and laser STAC analysis before and during exposure to hypotonic shock or saline with isotonic concentration.

Figure 12 presents the dynamic features of the procedure described above, measured by speckle dynamics drift velocity. Whereas for the specimen with control tonicity negligible drift was seen, speckle drift increments as shown in Fig. 12(B), demonstrates rapid process, with the swelling process onset at about 5 seconds post exposure to hypotonic solution. The drift observed in the above-described experiment was in the x direction due to the positioning of the camera with respect to scattering origin. A radially symmetric swelling process is expected to cause a corresponding radially symmetric speckle correlation drift, but the field of view of our camera could only collect a small segment of the whole speckle pattern. Nevertheless, considering the location of the camera with respect to scattering center axes, the two-dimensional picture of such radially symmetric speckle drift should occur along the x axis as seen in Fig. 12(A).

 

Fig. 12. Speckle correlation peak drift presented in 2D trace of both control (red) and hypotonic shock (blue). A. directionality of motion is displayed where drift is on the x axis due to the choice of camera location. B. dynamic behavior for the cases of control and osmotic shock.

Download Full Size | PDF

As mentioned above, image analysis tools could also be applied to trace the on-going swelling process. To further demonstrate the information that may be extracted from the STAC analysis method, one may focus on the differential increments of spheroid diameter either measured by image analysis or by speckle correlation drift measurement. A comparison of differential measurement is displayed in Fig. 13. The normalized process of incremental size growth as detected by image analysis is shown in Fig. 13(A) where the largest swelling step is shown to occur in the interval 0-10 seconds post exposure. During the following 60 seconds, the incremental swelling rate decreases with time.

 

Fig. 13. Differential size changes in spheroid post hypotonic shock, as measured by A) microscopic image analysis software and B) speckle STAC peak drift (with running average over 3 consecutive points)

Download Full Size | PDF

In Fig. 13(B), normalized speckle drift parameter trace is shown with what appears to be enhanced details of the swelling process. The additional information may be attributed to (1) the higher throughput capacity of the STAC method in respect to the slow rate of image acquisition in the present work and (2) the sensitivity of the STAC method to cell motion taking place inside the spheroid. To date, little is known about the exact structural dynamics taking place inside a spheroid undergoing osmotic shock. The results we obtained by applying the STAC method to the measured speckle pattern clip are not fully supported by any other current method. Further efforts should be made to validate the structural dynamics implied by the results presented in Fig. 13(B).

4. Discussion and conclusions

In this work, a single-scattering approach of speckle spatial-temporal auto-correlation (STAC) model was introduced as potentially exploitable for the analysis of three dimensional multiple scattering objects such as spheroids. Although limited to a small angle scattering regime, the STAC procedure has been found to yield encouragingly accurate predictions as to motional parameters inside multiple scattering spheroid models. The main feature of the STAC analysis is extraction of the directionality of inflation motion inside the spheroid, in contrast to non-directional velocity information that is usually extracted using DLS methods. Also, the level of order-disorder ratio in the motion may be deduced from the specific decay rates of features within spatial-temporal speckle auto-correlation distribution. Exploitation of this method on actual microscopic specimens was experimentally demonstrated by introducing a phantom which contained polystyrene beads as a three dimensional (3DS) scattering object. The internal motion induced in the phantom was numerically modeled to reproduce the speckle dynamics of the scattering pattern. Then, STAC analysis was applied on both experimental and simulated speckle patterns.

As mentioned above, the validity of the proposed procedure is expected to be limited to relatively small angles. The exact theoretical and experimental implications of these limitations are to be investigated in further study. Yet, there are two important aspects by which suggested STAC may add information not available by means of previous work [14,15]. The first, is the Synchronic Motion Parameter (SMP) presented in Eq. (29) which yields information regarding the synchronization of internal motion within the investigated object. The second is ${\alpha _{axis}}\; $ which concerns the directionality of this internal motion of the cells composing the 3DS. Both aspects could be tested by examining the experimental data obtained from the phantom 3DS. This data, which may be seen in Fig. 8, includes images of the internal structure of the 3DS object and of its corresponding speckle pattern as well, taken at a definite angular area. As depicted in Fig. 9, both measured and simulated speckle patterns where subjected to STAC analysis, thus applying a spatial-temporal correlation peak-tracking approach. The expected values for SPM and $\textrm{}{\alpha _{axis}}$ parameters, based on object plane images, where compared to the measured values as obtained by measured speckle records. As may be seen in Figs. 9(B) and 9(E), the ${\alpha _x}\textrm{}$ directionality parameter for the measured (Fig. 9(B)) and simulated (Fig. 9(E)) speckle patterns, closely resemble each other, as the slopes of the STAC distribution are almost the same. This holds true also for the comparison between experimental and synthetically predicted synchronization parameter SPM applied on the mentioned speckle data. As may be seen Figs. 9(C) and 9(F), the curves for DLS and STAC correlation functions deviate from each other for both experimental and synthetic speckle dynamics, which indicates an agreement between experimental results and numerically predicted parameter SPM, where the values of the measured $SP{M_{meas.}} \cong 0.68$ and numerically predicted $SP{M_{pred.}} \cong 0.63\; $ of the internal motion in the 3DS phantom are quite close to each other.

The above results obtained with a synthetic phantom 3DS, demonstrate the potential of the STAC approach as a method that could characterize both directionality and synchronization of internal motion within multicellular objects. In addition, the STAC method was applied on a living spheroid subjected to hypotonic shock while both STAC and image analysis, were compared. As was evident in STAC results, the size of the spheroid was elevated in a process that spanned dozens of seconds with temporal resolution of about a hundredth of a second. While both STAC and image analysis expressed similar dynamic behavior of the spheroid, the STAC analysis presented a significantly more detailed, fine structure that is yet to be understood and validated in future investigation. To conclude our discussion, we point out the potential of the proposed method as a fast, simple, and affordable method for the extraction of internal dynamics in 3DS.