Optical Tweezers with Integrated Multiplane Microscopy (OpTIMuM): a new tool for 3D microrheology

We introduce a novel 3D microrheology system that combines for the first time Optical Tweezers with Integrated Multiplane Microscopy (OpTIMuM). The system allows the 3D tracking of an optically trapped bead, with ~ 20 nm accuracy along the optical axis. This is achieved without the need for a high precision z-stage, separate calibration sample, nor a priori knowledge of either the bead size or the optical properties of the suspending medium. Instead, we have developed a simple yet effective in situ spatial calibration method using image sharpness and exploiting the fact we image at multiple planes simultaneously. These features make OpTIMuM an ideal system for microrheology measurements, and we corroborate the effectiveness of this novel microrheology tool by measuring the viscosity of water in three dimensions, simultaneously.

the position of the trapping beam waist and therefore of the optically trapped particle is also altered; (ii) when it is achieved by moving the sample holder, disruptive vibrations may overshadow the thermal fluctuations; and (iii) piezo scanners are commonly limited to speeds ≤ 100 Hz 14 , which would significantly hinder the highest accessible experimental frequency to which the viscoelastic properties of the material could be deduced. An alternative approach that doesn't suffer from any of these problems is to use a QPD. The Gouy phase shift obtained from a QPD can allow beads to be effectively tracked in 3D over small distances in z 15,16 . This has been combined with optical trapping to allow for 3D imaging of soft structures, to great effect 17 . However, QPD's can be challenging to calibrate, and only operate over a short range in z. In the literature, several other imaging techniques have been reported that also extend the tracking of a bead to 3D and, in principle, they could be employed for MOT; these include using electrically tunable lenses 18 , interferometry 19 , and stereo-microscopy 20,21 .
One other option is to use holographic microscopy, as demonstrated in the microrheology system of Cheong et al. 22 . This system gives excellent microrheology results, but in doing so replaces the standard illumination with a HeNe laser to produce the holographic image. This makes the system less convenient for probing ultra-local properties in biological systems, where it may be useful to simultaneously use transmission images to see how close the probe bead is to e.g. cells 6 .
The system we present here is based upon multiplane microscopy. Multiplane microscopy is a potent imaging technique that allows us to capture images at multiple different focal planes simultaneously. Conventionally, multiplane microscopy techniques can be classified into two categories: (i) those that split the image with beam splitters and then use lenses to change the focal depth for each beam path 23 , and (ii) those where diffractive optical elements are employed 24 . In this work we have used the latter approach to simultaneously image a single optically trapped microsphere at nine different focal planes and thus perform microrheology measurements in three dimensions. A significant advantage that OpTIMuM has over many of the aforementioned 3D imaging techniques is its simplicity, as it requires only the addition of passive optical components into the optical path between the microscope and the camera. Previous multiplane microscopy work by this group used a look-up table approach to track sub-diffraction limit sized beads in 3D 25 . This is a procedure that, in common with other look-up table approaches requires the generation of a calibration data set on a bead and sample system which must match the optical and the physical properties of the bead being used for measurements. These properties include, bead size, bead refractive index, fluid homogeneity and illumination. Generally the measurements will be carried out using the same tracer bead that the calibration data set was made with.
The ability to gain spatial information in the z-direction makes OpTIMuM a valuable tool for a variety of studies throughout the applied sciences. These include, but are not limited to, exploring the properties of 3D cell cultures 3,26 , elucidating the mechanical proprieties of materials with highly anisotropic structures 27 , and investigating the complex interactions at the solid-liquid interface of planes and/or biochemically treated surfaces 28,29 .
Our approach to 3D localisation may equally be applied to studies of complex interactions and processes involving cells, macromolecules and pathogens [30][31][32][33] , where the combination of high frame rates, no moving parts and compatibility with optical trapping would make it highly effective.

Methods
The system was designed using an Olympus IX73 inverted wide field microscope with an Olympus 60 × water immersion lens (Olympus UPLSAPO60XW) having a numerical aperture (NA) of 1.2. Images were taken in transmission, using a green LED (Thorlabs LED-C13) for the illumination and then filtering the images using a narrow bandpass filter (Thorlabs, FWHM = 3 nm) centred on 532 nm. The narrow filter is required to minimize chromatic aberrations due to the gratings. Images were captured using a Hamumatsu Orca Flash 4.0 camera at 67 frames per second. To test the capability of the z-localisation we used a piezo microscope stage (Mad. City Labs Nanodrive 85) to move the stage in z. We worked in transmission rather than fluorescence to avoid photobleaching the optically trapped particle, which would have had increased the complexity of the z-localisation process. A schematic diagram for the system is shown in Fig. 1a. For the optical trap, we used the output of a 1064 nm laser (Opus, Laser Quantum). In our system, as in most optical tweezer platforms, a beam expander consisting of two lenses (L1 (30 mm) and L2 (200 mm)) ensured that the trapping laser over-fills the back aperture of the objective and makes best use of its NA. By adjusting the position of lens 2 relative to lens 1 via a manually operated Thorlabs CT1translation mount, any bead trapped by the laser would move in the z-direction relative to the planes of the imaging system. Similar configurations have been employed in optical tweezers experiments previously 29 . To avoid the effects of laser heating on the sample we used a 1064 nm laser which is only very weakly absorbed by water and kept trap power ≤ 20 mW throughout.
The multiplane system used in this work is similar to the one described in refs. 25,34,35 . To briefly summarise, we have attached a 4f image relay system (consisting of two 300 mm lenses) to the camera output port. The 4f system allows us to place the multiplane grating in the telecentric position and thus maintain a consistent level of magnification in each of the imaging focal planes. The multiplane grating itself is a quadratically distorted diffraction grating etched into a quartz substrate by Photronics UK Ltd. A schematic representation of the etch pattern is shown in Fig. 1b, although the real gratings have a much shorter etch period than shown. The optics underpinning the operation of the multiplane system are described at length in ref 36 , but in simplest terms the object planes from different depths in z are spatially separated on a single camera. A single grating may be used to separate an image into three sub-images, corresponding to the m = 0, ± 1 diffraction orders. By using two gratings with orthogonal etch patterns, it is possible to split a single image into nine different sub-images, which due to the quadritic distortion of the gratings have each corresponding to a different image depth, as shown in Fig. 2a. In our system we have used a relay and grating combination that gives plane separation of Δz = 0.88 µm, a short explanation of why this spacing was chosen can be found in the "Z-localisation" section where we discuss how sharpness may be used to localize the bead in z. The plane spacing was measured by moving fixed beads of The image processing is represented schematically by the flowchart shown in Fig. 2b. To summarise, raw images captured by the camera are segmented into the nine different sub-images each corresponding to a z plane = n Δz, for n = ± 4,3,2,1,0, where the central plane is n = 0. Each segmented image stack is then processed with the "multithresh" function in Matlab [Matlab 2019b; MathWorks, Natick, MA]. This uses Otsu's method to find a threshold value which maximises inter-class intensity variance and minimises intra-class intensity variance 37 . This threshold is then applied to the segmented images to create a black background, and the x, y coordinates calculated by performing centre of mass analysis. The segmented images without the threshold are retained, a background level is subtracted, and the absolute value of the pixel intensities taken. These are then used to calculate the z-position of the particle as described in the "Results" section.
The underlying principles of the microrheology analysis performed on the particle trajectory is extensively described in references 9 and 10 . Briefly, we calculate the normalized mean squared displacement (NMSD) of the bead in x, y and z directions, which is related to the normalized position autocorrelation function (NPAF) by the following equation: where κ = trap stiffness, defined as κ x = k B T/ < x 2 > (where x may be replaced by y or z, T is the absolute temperature and k B is Boltzmann constant) and i = x, y, z . Bead radius is given by r, τ is lag-time and η is the unknown fluid viscosity that could be measured as described in the "Results" section. Note that the last equality of the above equation is valid only in the case of Newtonian fluids.
To test our z-localisation method, we used a bead embedded in a stiff agarose gel such that its thermal motion in z is negligible compared to the 50 nm steps of the piezo stage. This was achieved by dispersing a low melting point agarose (TopVision, R0801) in warm (> 65 °C) distilled water and then adding a very dilute suspension of beads into the melt, before pipetting a small amount of sample onto a microscope slide and letting it set at room temperature. The sample was then mounted on the microscope, a single bead brought into focus, and measurements performed by moving the sample through a known distance via the piezo z-stage. Several different batches of beads were used with a range of different diameters; d = 2.1 µm (Polysciences Inc, 19508), d = 3 µm (Polysciences Inc, 17134), d = 4.5 µm (Polysciences Inc, 17135), nominal d = 5 µm (ThermoScientific 35-2, these beads were actually highly polydisperse, d = 5-8 µm ) and d = 9 µm (Invitrogen, N37464).

Results and discussion
Z-localisation. In order to perform microrheology measurements, an accurate localisation of the trapped bead is needed. In x and y directions, this is typically achieved by thresholding the image and then taking the centre of mass of the bead intensity profile in each dimension, as defined by Eq. (1). Here I xy is the intensity of the pixel with indices x and y in each dimension and the sum is over all pixel indices.
However, as shown in Fig. 3a, for a colloidal bead imaged in transmission, the sum of the image intensity has poor signal to noise ratio (blue curve) and does not peak when the bead is sharply in focus. Therefore, the centre of mass method used in x and y cannot simply be extended to the z-dimension, and a different approach must be used, as described below.
(1)  (2), which is a single metric that peaks when there are minimal aberrations and falls off rapidly as aberrations increase 38 . Our system has narrow spectral filtering, uses high quality optics and is carefully aligned to minimise aberrations, so the shape of the sharpness curve will be dominated by defocus. In Fig. 3a we show the sharpness curve for a bead being scanned in z (orange curve). Steps on the blue background relate to x and y localisation, whereas those on the gold background relate to z localisation; those on the green background represent initial processing steps common to x, y and z.  34 , a fixed bead was repeatedly moved by a known distance in z and the sharpness in each plane measured. This created a matrix of known sharpness values for each image plane in each stage position that was then used as a calibration dataset. By using statistical techniques based around maximum likelihood estimation, it was possible to calculate the position of this bead by comparing the measured Sharpness www.nature.com/scientificreports/ values to the calibration set and finding the z value most likely to produce those Sharpness values. This offered excellent precision, when used for sub-diffraction-limit sized beads rather than the 1 -10 µm diameter beads commonly used in optical trapping based microrheology. Similar "look up table" based approaches have also been used for larger beads, and based on other metrics such as radial profile rather than sharpness [39][40][41] . However, regardless of the metric used, such approaches strongly rely on the ability of the user to generate their calibration by simply scanning their objective through z while holding the bead stationary. This is achievable for magnetic tweezers measurements, where the position of the trapped object can be controlled independently of the imaging plane, but for the majority of optical tweezers systems imaging and trapping are performed via the same objective lens which makes this far more challenging as it is not straight forward to move the bead relative to the image plane. A possible alternative would be to generate a calibration set using a different bead from the one which will be used for measurements. However, this depends on the assumption that the behaviour of the metric in the calibration set will be identical to those of the measurement. In reality, optical properties vary greatly depending on the depth being imaged at and the inhomogeneity of the sample microstructure (this is especially true of complex biological samples), and there will always be a degree of polydispersity in bead size (the Sharpness function is strongly dependent on bead size, as shown in Fig. S1). Therefore, instead of using Sharpness directly as a metric for a look-up-table, we have explored the use of Sharpness as a weighting function for a centre of mass calculation in the z direction (similarly to how Intensity is used in the x-y plane) returning Eq. (3), where z C.S is the centre of sharpness.
In Fig. 3c we show a comparison between z C.S (red line) and the actual z position taken from the piezo controller (black dashed line), for a ~ 3 μm diameter bead fixed in agar and imaged in transmission with 15 ms exposure, as the stage is moved in the z direction in 50 nm increments. It is clear that z C.S significantly underestimates the degree of motion of the bead, as the gradient of z C.S (~ 0.2) is much lower than the gradient of 1 exhibited by z stage . This is because with only nine planes in z we are significantly undersampling the function (see Fig. 3b), and hence biasing z C.S to be too close to zero.
An alternative approach is to identify the plane that has the highest Sharpness value at any given time (z sharpest ) and take the z value of this object plane relative to the mid-plane to be the bead position. This method is accurate, but offers very low precision, equal only to the plane spacing. This can be improved by finding the plane with the second highest Sharpness value (z 2nd-sharpest ), and assuming that the bead must lay between these two planes, but closer to the sharpest plane. We define this as the 'Sharpest Plane' approach and z S.P as in Eq. (4). Equation (4) returns the same values as z sharpest but shifted by ± Δz/4 depending upon whether z 2nd-Sharpest is the plane above or below z sharpest This results in even step sizes of Δz/2 as the bead moves in z. Although z S.P is more precise than z sharpest , it still only returns a resolution of half our plane spacing, as shown in Fig. 3b by the blue step-wise line. This approach lacks the spatial resolution of < 100 nm essential for microrheology studies.
In order to achieve the required spatial resolution, we have combined the two methods outlined above. Specifically, we have used z S.P as a calibration to linearly rescale z C.S . We have taken some additional steps to avoid the effects of non-linearity when the bead has moved far from equilibrium, and to reduce biasing (which we detail in Sect. 3 of the SI), to perform a linear fit of the stepped z S.P data as a function of z C.S . For the data shown in Fig. 3 this returns a gradient of 5.26, which we take as a rescaling factor specific to a bead of this size under this illumination, i.e. z R.C.S = 5.26 z C.S . It can be seen that z R.C.S (green line) follows the actual stage position (black dashed line) very closely, as quantified by the residual plot. Notice that Fig. 3b was generated by moving the stage a total distance of only 4 µm. The calculated z position matches the stage position very well over this range, with a mean residual of ~ 40 nm. However, it must be highlighted that in actual hybrid optical trapping microrheology experiments where trap strengths are generally between 1 × 10 -6 and 1 × 10 -8 Nm −1 , an optically trapped bead will typically move less than ± 1 µm in any dimension during a measurement performed at room temperature. This is confirmed in Fig. 4a,b, where the scatter plots of an optically trapped bead of ~ 4 µm radius are reported. Notice that the bead travels only ~ 1 µm in any direction from the centre of the optical trap. Therefore, when considering the accuracy of this method as it will be used in microrheology measurements, it is more relevant to consider the error in detecting the bead position for displacements of + /-1 µm from the centre. In this region we obtain a mean absolute residual between z R.C.S and z stage of z stage −z R.C.S ~ 30 nm for the d = 3 µm bead shown in Fig. 3d.
In order for this method to be effective, the plane separation must be large enough to adequately sample the peak and drop off of the sharpness function. This is demonstrated in Figs. S2-S5 where we show plots of the same type as Fig. 3c,d, but for beads of different sizes. For all the bead sizes explored, there is a region in which z RCS scales linearly with stage position. For larger beads, within this central region the method is more accurate, and residuals are smaller. However, it must be noted that as bead size becomes larger, the range of z stage values over which z CS behaves linearly becomes narrower, even as the residuals within these regions decrease (see Fig. S6). This clearly demonstrates that one must optimize bead size for plane separation as larger beads will give less random noise but may introduce systematic error if a very weak optical trap is used, such that the bead regularly travels over a range wider than this linear region. For the Δz = 0.88 µm plane spacing employed in this work, we conclude that using a 9 µm bead may risk introducing systematic underestimation of z motion into MOT experiments, whereas the other beads tested would safely give linear behavior over the full extent of the optical trap. The ThermoScientific beads with diameters ranging from 5-8 µm offered a precision of  Fig. S4) within the ± 1 µm region, therefore, we have used beads from this batch for the measurements reported hereafter. The proposed methodology depends on moving the bead up and down at the start of the measurement to generate the z S.P function required to perform the in situ self-calibration procedure. This can be achieved in different ways depending on the nature of the sample under investigation. For instance, in the case of a bead embedded in a very stiff media capable of holding it in place (such as a gel), the use of any standard z-stage with step precision ~ Δz/2 or capable of continuous motion in z will allow this task to be accomplished. This is more complex for microrheology measurements of fluids, because the trapping laser and transmitted illumination light travel through the same objective. Therefore, moving the sample or the objective would not change the relative position of the trap (i.e., the bead) to that of the imaging planes. However, by moving lens L2 (as described in the method section and indicated in Fig. 1a) it is possible to adjust the position of the laser beam waist independent of the imaging plane. This allows us to move the bead over a distance of ± 4Δz, which is more than sufficient to obtain z C.S values at a range of different values of z S.P , thus allowing a correct rescaling of z C.S to generate z R.C.S . This is demonstrated in Fig. S6, which shows z C.S , z S.P and z R.C.S for the self-calibration step followed by the initial part of a microrheology measurement.
3D Microrheology with optical trapping. In order to corroborate our method, we performed preliminary microrheology measurements in water, a Newtonian fluid. Here we continue to operate at 67 Hz, which is adequate to analyse the viscous behavior of water. However, it important to highlight that, as there are no moving parts in the multiplane system, this speed is only limited by the camera acquisition rate. For instance, by using gratings with a larger etch period, or a shorter relay with gratings with higher curvature, we could reduce the field of view and thus the number of pixels in our images. We could then operate at ~ 500 Hz using the same camera. Figure 4a shows a 3D scatter plot of the trajectory of a bead confined within a 3D optical potential. The trajectory has been drawn by analysing 100,000 individual frames each recording 9 simultaneous z planes imaged during the measurement. This allows the visualization of the entire volume of the optical trap, which has a prolate spheroid shape stretched along the z-axis, due to the inherently weaker trap strength in the axial direction 43,44 . The spatial distribution along the three axial directions has a full-width-half-maximum of 0.45 µm, 0.45 µm and 0.76 µm in, x, y and z, respectively (Fig. 4b). The x, y and z coordinates versus time are plotted in Fig. 4c-e and do not show any significant drift, indicating a thermally and optically stable system. This is confirmed by looking at the Allan deviation of the bead position in each dimension which also shows no drift on the timescale of the experiment (see Fig. S7). Calculating the trap strength gives us κ x ∼ = κ y ∼ = 1 × 10 −7 Nm -1 and κ z ∼ = 4 × 10 −8 Nm -1 . The particle trajectory can be further exploited to evaluate its normalised mean www.nature.com/scientificreports/ square displacement (NMSD) 10 as a function of elapsed time (τ) for each spatial dimension (x, y, z), as shown in Fig. 4f. We can see that the NMSD curves evaluated in the x and y directions overlay almost perfectly on top of each other, as expected for a spherical particle in an isotropic fluid and a well aligned optical system. The NMSD curve in the z direction, however, is offset in the time axis because of the weaker trap strength. This offset can be compensated for by following the analytical approach introduced by Tassieri et al. 42 , and plotting either the NMSD or the NPAF against a dimensionless lag-time, τ* = τκ/(6πrη s ), which takes into account the bulk viscosity of the solvent (η s ), as shown in Fig. 4g. We see that all the curves collapse onto a single master curve, as we are probing the same physical phenomenon, i.e. the Brownian motion of the same bead suspended in the same isotropic Newtonian fluid. Notably, this representation allows a direct measurement ('at a glance' 42 ) of the fluid relative viscosity ( η r = η/η s ), defined as the ratio of the fluid viscosity to that of the solvent. The relative viscosity will be the τ* value at which NPAF(τ*) = e −1 ; the latter is marked by the solid horizontal line in Fig. 4g.
Notably, the abscissa of the intercept of this line with all three NPAFs occurs at τ* ≅ 1, meaning that the measured relative viscosity is 1 as expected for pure water. Measurements were repeated on 5 beads with d ~ 6-8 µm and obtained η x r = 1.02 ± 0.04, η y r = 1.01 ± 0.05, η z r = 0.98 ± 0.07. This corroborates the overall effectiveness of the proposed experimental procedure.

Conclusion
In this work we have introduced a novel 3D microrheology system named "OpTIMuM" that combines, for the first time, optical tweezers and multiplane microscopy. OpTIMuM allows us to perform particle tracking microrheology in three dimensions with ~ 20 nm accuracy along the laser beam axis, without the need to generate a 'look up table' and without knowledge of the bead size, the optical properties of the suspending medium, nor the use of any high precision positioner during the measurement. A straightforward manual adjustment to the position of a lens prior to measurement is all that is required to calibrate the system. Notably, since multiple planes are acquired simultaneously and not through a scanning process, the acquisition rate of our system is limited only by the camera frame rate. Using this system we were able to obtain the viscosity of water in all three spatial dimensions simultaneously. OpTiMuM provides a fast and effective method for exploring the full 3D micro-environment of an optically trapped bead. We envisage it could also be useful when measuring the forces involved in cell-surface interactions studies or when looking at force extension relationships of single molecules, and could open new routes to topographical, chemical and bio-mimetical surface engineering studies for biomedical applications.