Using MATLAB to Achieve Nanoscale Optical Sectioning in the Vicinity of Metamaterial Substrates by Simulating Emitter-Substrate Interactions

The imaging performance of a conventional far-field microscope is restricted by the so-called diffraction limit Born &Wolf (1997). For an imaging wavelength of λ this corresponds to a maximum attainable resolution of approximately λ/4 (laterally) and λ/2 (axially) when using a high numerical aperture objective lens. With recent advances in molecular biology some of the most interesting questions require a resolution beyond this. Higher resolution has been typically achieved using lower wavelength probing/scanning fields, or resorting to near-field scanning techniques (e.g. Scanning Near-field Optical Microscopy (SNOM) Betzig et al. (1991)). Unfortunately these come with their limitations and undesirable side effects, which prove to be limiting when it comes to studying molecular dynamics in their native environment. Furthermore, the use of fluorescent labels for studying specific processes and background suppression is desirable in many applications. Far-field optical microscopy still remains the only practical non-invasive technique suitable for imaging fast dynamics over extended periods and distances of interest with minimal perturbations to the system. Fortunately, far-field optical microscopy capable of imaging beyond the conventional diffraction limit for life science applications has seen many advances over the last few decades. Our current arsenal of techniques include commercially well developed techniques such as confocal laser scanning microscopy Cremer & Cremer (1978) and Total Internal Reflection Fluorescence (TIRF) microscopy Axelrod (1981)1. The improved resolution with these techniques is in itself however only by a factor∼ 2 and often not sufficient inmany cases. More powerful techniques such as Fluorescence PhotoActivatable Localization Microscopy (FPALM), Stimulated Emission Depletion (STED) Microscopy, Structured Illumination


Introduction
The imaging performance of a conventional far-field microscope is restricted by the so-called diffraction limit Born & Wolf (1997). For an imaging wavelength of λ this corresponds to a maximum attainable resolution of approximately λ/4 (laterally) and λ/2 (axially) when using a high numerical aperture objective lens. With recent advances in molecular biology some of the most interesting questions require a resolution beyond this. Higher resolution has been typically achieved using lower wavelength probing/scanning fields, or resorting to near-field scanning techniques (e.g. Scanning Near-field Optical Microscopy (SNOM) Betzig et al. (1991)). Unfortunately these come with their limitations and undesirable side effects, which prove to be limiting when it comes to studying molecular dynamics in their native environment. Furthermore, the use of fluorescent labels for studying specific processes and background suppression is desirable in many applications. Far-field optical microscopy still remains the only practical non-invasive technique suitable for imaging fast dynamics over extended periods and distances of interest with minimal perturbations to the system. Fortunately, far-field optical microscopy capable of imaging beyond the conventional diffraction limit for life science applications has seen many advances over the last few decades. Our current arsenal of techniques include commercially well developed techniques such as confocal laser scanning microscopy Cremer & Cremer (1978) and Total Internal Reflection Fluorescence (TIRF) microscopy Axelrod (1981) 1 . The improved resolution with these techniques is in itself however only by a factor ∼ 2 and often not sufficient in many cases. More powerful techniques such as Fluorescence PhotoActivatable Localization Microscopy (FPALM), Stimulated Emission Depletion (STED) Microscopy, Structured Illumination Microscopy (SIM), and variations thereof -see e.g. review Hell (2007) -have been developed which can offer lateral resolution improvements by factors of up to ≈ 20. Despite the great success of these techniques they struggle with achieving sufficient temporal resolution in many cases. They also struggle with achieving axial super-resolution, which when achievable usually requires elaborate modifications to the setup. Currently there exists to our knowledge no far-field optical microscopy technique that is capable of studying axial dynamics with < 20nm resolution in live cells without either significantly perturbing the sample or requiring an elaborate and/or costly microscopy setup. Here we outline the computational aspects associated with a technique -Metamaterial Substrate Modified Fluorescence (MeSuMo) Microscopy -we have recently developed that allows for dynamic imaging with axial superresolution which is compatible with standard epifluorescence microscope setups. As with majority of optical imaging techniques currently used in the lifesciences, MeSuMo is essentially a fluorescence microscopy technique, requiring the objects of interest to be labelled with suitable fluorescent emitters. Compared to other techniques capable of achieving comparable resolution (e.g. Kanchanawong (2010)), MeSuMo microscopy is, once optimized, experimentally very simple to implement and suitable for fast dynamic studies on live cells. It relies on using microscope slides coated with an optimized metal-dielectric coatings (metamaterials) which modify the emission spectrum of fluorophores in a manner that is dependent on their separation from the surface of the substrate. A spectral analysis over a sampled area allows one to infer the average separation of the emitters in the studied region. In this chapter we will not aim to provide an exhaustive description of the underlying principles and details associated with the technique, but rather emphasis how the data analysis can be performed in a MATLAB environment. We note that there exists an extensive underlying theoretical body of work in the various areas of classical and quantum physics for modeling the near-field electromagnetic interaction with plasmonic structures and the photophysics of fluorophores, which the interested reader may learn more about through the referenced literature. To maintain the focus of the chapter we will thus in most cases only present qualitative interpretations of the effects when relevant and quote the essential results (equations and integrals) in a form that can be readily entered and computed using MATLAB. The chapter is divided into five sections that include: (1) modeling the metamaterial coated substrates, (2) fitting the experimental data to the model, (3) modifying the model to obtain best fit results, and (4) an example of the technique being implemented. We conclude the chapter in section (5) with a discussion of future additions to the technique we are actively working on.

Modeling the metamaterial coated substrate
An optical metamaterial is an artificial material typically consisting of subwavelength scale metal and dielectric components. In the most general context its key feature is, in a nutshell, that it has a unique response to certain electromagnetic fields that can not be found in naturally occurring materials. This typically consists of negative refractive properties -namely that an incident field would refract in the opposite direction from the normal-axis as compared to a conventional naturally occurring material. Numerous interesting and potentially useful effects may result such as strong enhancements of the Purcell factor for nearby emitters Jacob et al. (2010), which allow for the realizations of devices such as SPASERs Bergman and Stockman (2003) and lasing SPASERs Zheludev et al. (2008). Here we focus on a relatively simple case of layered metallic-dielectric layers each having a thickness much smaller than or comparable to the wavelength. Depending on the precise thicknesses of the layers as well as their individual response properties such structures may exhibit Fabry-Perot resonances, resonant photon tunneling and other unusual transmission/reflection properties Belov and Hao (2006); Darmanyan and Zayats (2003); ; Ramakrishna et al. (2003). It turns out that many of these properties can quite intuitively be modeled analytically using MATLAB due to the matrix formalism often used to describe them. In the context of the technique we will discuss there are two properties of interest to us. Firstly, there is the complex reflection coefficients of the structures. These determine the local field at the site of the emitter and hence also the excitation and decay rates. Secondly there is the dispersion relation of the supported modes of the structure which will give us physical insight into the type of excitations the emitters couple to at different distances and frequencies, and thereby allow us to tune our structures to obtain optimal results. Of particular relevance for the latter will be a so-called "cut-off energy" which exists in asymmetric structures at a given energy, where there is a transition from a bound to an unbound SPP mode Burke and Stegeman (1986). For the case of the transmission/reflection coefficients one simple way to model the structures is to use transfer matrices which describe the propagation and attenuation in each layer Born & Wolf (1997). As mentioned such calculations are particularly well suited for MATLAB where, once the matrices are defined, the coefficients for arbitrary combinations of layers and structures can readily be determined and compared. The transfer matrices are generally defined in the Fourier domain parallel to the surface of the layers (k x ,k y ) and the real domain normal to the layers (z). This also allows for easy analysis of the contribution from different evanescent field components which are dominant in the near-field (distances smaller than about the wavelength). Once one has written matrices for the transmission through a given interface (t i,j ,wherethesuperscriptsi and j denote the regions on either side of the interface) and the attenuation through a certain thickness of a given material (p i )i nt e r m so ft h e transverse wavevector(s) and -for the latter case -the thickness of the layer, one can obtain the total transmission function (or optical transfer function) of an n-layered structure t n by t n = Π n−1 i p i t i,i+1 p n .( 1 ) Details on the implementation of transfer matrices can be found in most standard classical optics texts Born & Wolf (1997) and will not be elaborated on here. We only mention that the formalism will, for not too large transverse wavevectors, also be applicable for subwavelength laterally structured layers (e.g. layers containing split-ring, horse-shoe shaped, or fractal resonators). In such cases one would need to define suitable effective anisotropic permittivity and permeability matrices for them, e.g. in a 2D system one may have a permittivity epsilon(a,b) and permeability mu(a,b) where a and b are the x and z components required to describe the effective response properties of the structure. The cases of two polarizations (perpendicular magnetic and electric fields) can be treated separately by choosing the suitable Fresnel coefficients for the matrix elements of t i,i+1 . Alternatively one can also use a "brute force" approach for not too complicated structures and include the complete equations (e.g. as a function) into MATLAB. For example, for the case of a four layered structure (layers indexed by i, j, k, l)thiswouldbe: and where ε i , μ i and d i are the permittivity, permeability and thickness of the i th layer. The total wavevector in the ith layer is given by where ω is the frequency of the electromagnetic field and c is the speed of light in vacuum. Also a normalized in plane wavevector is defined as u = k x /k 1 ,wherek x is the in plane wavevector component which is independent of the layer. The value of ε for the metal will in general be complex and frequency dependent, and is most easily obtained by fitting a high order polynomial to experimentally measured data (e.g. Drachev (2008) for silver). Alternatively for not to thin metal layers one may to a good approximation use a Drude-Sommerfeld model (e.g. Kittel (1995)) where the parameters can later be systematically varied to account for e.g. changes in temperature Elsayad and Heinze (2010b).
To determine the dispersion relation and the cut-off energy of a given structure one may also make use of the natural matrix calculation approaches of MATLAB. To do so one firstly writes a general expression for the parallel electric fields in each layer. For a layer-i extending ,w h e r ek zi is the out of plane wavevector in the layer and A are at this point arbitrary complex coefficients. One may then write the corresponding magnetic induction fields for each layer using Maxwell's equations (i.e. by taking the normal derivative), so that for each layer one has two equations with two unknowns. This is done for each layer i = 1, 2...N. Finally one constructs a 2N × 2N with the coefficients for each layer constituting the entries in each row. The determinant of this matrix will then give the dispersion relation ω(k x ) -s e e e.g. Dionne et al. (2008). Since k x is complex one can not simply assign a solver such as fzero to the task. It is it turns out often necessary to evaluate this by brute force -i.e. letting MATLAB evaluate the determinant for a given complex transverse wavevector and checking how close it is to zero. In general it is important to provide a good starting guess of Re(k x ) and Im(k x ) for each frequency. For the cases of interest, where the frequencies of interest are close to the cut-off frequency, a good starting guess for the former is the light-line Re(k x )=nω/c of the medium with refractive index n in which the mode becomes unbound, and then to progressively search further into the evanescent region, i.e. Re(k x )=nω/c + δ with δ ∼ 0.01nω/c. For the imaginary part Im(k x ) -which diverges at the cut-off -it is most efficient to employ an algorithm which calculates the gradient of the determinant between two nearby values of Im(k x ) [i.e. taking the difference in the determinant value and dividing by the difference in Im(k x )], and subsequently search for larger/smaller Im(k x ) values if the determinant is smaller/larger than zero. The stability of the solution should be checked to assure Lim δ→0 {ω(k x + δ) − ω(k x − δ) → 0}. The above mentioned approach typically requires that Im(k x ) is well behaved and monotonous in the vicinity of the cut-off, which is fortunately often the case. Since the equation has to be satisfied for both Re(k x ) and Im(k x ) optimizations of both of these have to be performed simultaneously. For completeness an analysis can also be performed for transverse electric (TE) as well as transverse magnetic (TM) polarizations, although for the cases of the very thin non-magnetic plasmonic structures only the latter are usually relevant. For the case of highly asymmetric structures of few layers, coupling between the interfaces may be negligible and a single plasmonic mode at one interface can to a good approximation be assumed solely responsible for the dispersion of the cut-off mode Burke and Stegeman (1986). In this case the cut-off energy can for a non-magnetic structure be estimated by E c =hcε j ε k [ε i (ε j + ε k )] −1 ,whereε i is the permittivity of the medium where the cut-off occurs, and ε j and ε k are the permittivities on either side of the interface at which the mode is originally localized. The principle of MeSuMo requires that the cut-off of the structure falls within the emission spectrum of the emitter. Since many common fluorophores and dyes have a fairly broad emission spectrum this criteria can usually be met quite easily and is generally quite robust. To calculate the change in intensity at a given frequency one firstly defines the change in the decay rate of a fluorophore as a function of distance and frequency. This may (at not too small distances) be modeled quite well in the dipole approximation using the classical Chance-Prock-Silbey (CPS) model Chance (Prock and Silbey). Under these conditions the decay rate in the +ẑ and −ẑ direction (away from and towards the substrate respectively) for perpendicular (⊥) and parallel (||) electric dipole like emitters can be obtained by: in which the integrands are where q is the quantum yield of the emitter, the superscripts P(S) denote the transverse electric(magnetic) cases, and the z and ω dependence ofΓ The integrals can in most cases quite easily be evaluated numerically (e.g. using quad). It may be the case that the structures have resonances (as can be determined by plotting the integrands as a function of u), in which case one may want to re-write the integrals to assure one integrates in a complex plane around them and account for the pole(s) in the conventional fashion. Integration can in almost all cases safely be cut-off at u ∼ 1000, with the contribution beyond u ≈ 100 becoming negligibly small for z > 10nm. In what follows we will sometimes make use of the notation ] ω ij to signify the isotropic average total decay rate from state i to j separated in energy by Besides the modification in the decay rate which has been described above, the measured emission intensity of an emitter will also depend on it's excitation rate. In particular the measurable emission intensity for the transition from the excited state a to a lower state b, which corresponds to a wavelength λ ab = 2πc/ω ab will be given by: where E ex (ω ab ) is the excitation field at the frequency ω ab ,andμ a is the dipole moment of the excited state a. Γ R ab and Γ a are the radiative decay rate for the transition a → b and the total decay rate of the state a. f ab is the Franck Condon coefficient between a and b, which is related to the rotational and vibrational configurations of the two energy states relative to each other, and can be thought of as the intrinsic favourability of the transition. For a realistic emitter there will in general be many other significant "down transitions" originating from the same excited state. For the decay of the state a to any other state x, one may thus write: Subsequent decay of all states x (and state b) to a common ground state is assumed to proceed non-radiatively and at a rate several orders of magnitude larger than the radiative decay rate. For the case of spontaneous emission we consider it can for all intense and purposes be assumed to be instantaneous. If we neglect direct resonant energy transfer between emitters (i.e. assume sufficient dilution), then the total excitation field with a frequency ω ab at the location of the fluorophore will be the sum of the external incident field E 0 (ω ab ),thereflection of this field from the structure re −2(1−u 2 ) 1/2 k 1 z E 0 (ω ab ) where r is the reflection coefficient, and the field generated from excitations in the material by all the fluorophores E r (ω ab , z),i.e.
The polarization for the reflected field can in most cases be taken as that of the incident excitation field. For the case of a single fluorophore E r (ω ab , z) is given by the integrands of equations (6) and/or (7), which are identical to the reflected field for a given transverse wavevector component u, and will henceforth be written as E 0 r (ω ab , u, z). To calculate the reflected field from all of the fluorophores at various distances from the substrate one can proceed in several ways. It is of no practical relevance in our applications if one requires information on the position and orientation of each and every contributing fluorophore. We thus can proceed only via an effective medium approximation. We will assume a density ρ of fluorophores which may only vary in the axial (z)direction.Inthecaseforρ(0 < z < z max )= ρ and ρ(z > z max )=0 one obtains by integrating over z ={0, ∞}: Alternatively one may represent the contribution from a collection of fluorophores by assigning a negative imaginary response function to the top layer (layer-1) in which they are immersed -i.e. set Im(ε 1 )< 0 for non-magnetic structures. Whilst this may be computationally simpler for constant densities, for the case of a non-trivial z dependence the calculations become significantly more elaborate. We thus proceed with the former method. For imaginary (1 − u 2 ) 1/2 k 1 the coupling fields are propagating and equation (12) can be reduced to the result for e.g. the field in a large semi-transparent 1D box filled with emitters. For real (1 − u 2 ) 1/2 k 1 the field will on the other hand decay rapidly with increasing z as one would expect. The contribution to I(ω ab ) from equation (12) will in itself however only be small even for large ρ (> 0.01nm −1 ). There will be an additional second order contribution to the excitation field and hence the emission intensity at ω ab from excited surface modes of energy ω ax where x = b.T h i s contribution is only significant because the SPP modes have very short decay times -typically Γ ∼O(10-100fs), such that their uncertainty broadening is relatively large -O(1 − 10eV). For ω ax ∼ ω ab , or more specifically |ω ax − ω ab |∼O (meV), they will contribute significantly to E R (ω ab , z). It follows that when the energy ω ab corresponds to mainly real values of (1 − u 2 ) 1/2 k 1 (evanescent fields) whereas ω ax to mainly imaginary values of (1 − u 2 ) 1/2 k 1 (propagating fields) this contribution becomes large near the surface as fluorophores far away from the substrate can effectively pump those nearby. To the lowest order the contribution can be written as: where and Γ ab (u, z) and E r (ω, u, z) are the corresponding values of the total decay rate and reflected field as a function of transverse wavevector [i.e. the integrands of equations (2) and (2) &/or (3)]. We note that equation (13) accounts for the energy overlap between a state of energy ω b and a continuum of states whose energy is given by the dummy variable ω β .T h e ω α integration is performed over a range of frequencies in the vicinity of the cut-off (typically corresponding to a wavelength range of 30-80nm), over which the energy range coincides with the emission spectrum of the fluorophore (i.e. are allowed transitions for the particular fluorophore). The weighting (Franck Condon) factors f αβ are empirically determined from the emission spectrum of an isolated fluorophore scaled so that the value at ω ab it is unity. We note that it is usually practical or necessary to convert this integration into a discrete sum with an energy spacing of ∆ω ∼ 1meV. Information particular to the type of fluorophore (e.g. energy level positions/separations and corresponding Franck Condon coefficients) may be implemented when constructing the summation. Since we will not focus on single molecule studies we assume an isotropic average of dipole orientations for all cases. This can most easily be done by setting E iso and/or Γ iso =( 1/3)Γ ⊥ +(2/3)Γ || as appropriate. We also note to calculate the final measurable emission intensity the integral for Γ R which appears in equations 9 & 10 (and is given by equations 6 & 7) should only be performed up to √ ε i ·NA, where NA is the numerical aperture of the measuring objective. It is possible to now predict the change in the measurable emission intensity for arbitrary z and fluorophores, near a metamaterial structure.

Fitting the experimental data to the model
Data can be acquired using a standard scanning epifluorescence microscope with a spectrometer and photomultiplier tubes. It is generally most practical to save images as .lsm files which can be directly read into MATLAB [In our setup each file contains a scanned image of the sample at ten different wavelengths (e.g. bins centered at λ = 490nm, 500nm, 510nm, etc.)]. For dynamic studies images at different time intervals can also be stored in a single .lsm file. Many spectrometers come with MATLAB drivers which allow for acquisition to also be controlled easily through MATLAB. Else for real time (video type) imaging it is necessary to save the data incrementally into a directory accessible to MATLAB. The spectral profile of the fluorophores of interest on an uncoated substrate immersed in a similar refractive index medium as the actual samples should also be available. The area under this spectrum should be normalized and extrapolated to a discrete matrix [e.g. fitted to a high order polynomial and evaluated at the discrete wavelengths indexed by i], and saved as a 1 × 10 matrix Ilambda0(i). Data analysis is performed on a pixel by pixel basis. For each pixel the spectrum area is normalized, fitted with a high order polynomial, evaluated at i discrete wavelengths, and written into a 1 × 10 matrix data(i). The ratio enhanobserved(i) = data(i)/Ilambda0(i) is then calculated. A 1 × j matrix zn(j) is defined which will represent the fraction of the emitters in the studied area at a distance corresponding to the index j from the interface. The task is now to solve enhan * zn-enhanobserved=0,which is in principle straight forward using a function such as e.g. linsolve. The result is best presented in a plot of zn as a function of z for each or a collection of pixels, which will show the distribution of fluorophore distances from the substrate at this particular lateral position. Depending on the sample studied and the information of interest this may be fitted with a suitable distribution function to obtain an average and spread in distances of emitters from the substrate within the analyzed region.

Optimizing the model and code
In many cases the initial fluorophore-concentration, metamaterial or coupling parameters may not be the best or correct parameters for the fit. This is usually evident when the results show an unexpected distribution for zn. Distances less than 10nm should be ignored as the approximations of the models (dipole approximation and local dielectric functions) are no longer valid at such small distances. The most significant source of uncertainty is the predefined and assumed constant value of the fluorophore concentration rho. To account for this one may redefine rho as being proportional to the integrated area of zn over z,a n d repeat the calculation of zn. Also to account for an axial variation in rho -as would be relevant for samples where the fluorophore density varies strongly with distance z -one may define rho(j) = rho * (a + zn(j)/norm(zn)),w h e r ea is normalization constant chosen so that on summation over the distances (i) one recovers the total intensity. In this case all subsequent calculations have to explicitly account for the additional distance dependence in rho (j), which significantly increases the computational time. This may be efficiently achieved by defining rho(j) in a smaller parameter space indexed by n that corresponds to dividing rho into components where rho(j) is constant and non-zero only over a specified region. Once zn is calculated in this manner the optimization may be repeated for the new rho until the the calculation of the distribution of zn is stable and smooth or varies on the expected scale.

Example
Whilst elaborate custom substrates may be designed with effective parameters which yield the desired cut-off energy and dispersion relation, it is often possible to use combinations of no more than 3-4 layers with no lateral structuring to achieve satisfactory results. Here we give an example where we use a sample consisting of a thin silver film 2 (15nm) deposited on a thick quartz substrate (ε = 2.13), which is subsequently coated with a high permittivity dielectric film (ε = 4.15). The sample is immersed in a medium with a refractive index larger than that of the quartz (ε = 2.4). The cut-off energy can be shown to occur within the emission spectrum of Alexa488 -a standard fluorescent marker with an emission spectrum peaking at a wavelength of around 528nm and extending all the way up to >550nm. To demonstrate the accuracy of the technique, fibroblast cells were grown on the laminin (thickness < 2nm, ε ≈ 2.15) coated substrate and a particular protein found at adhesion sites close to the substrate called Paxillin Schaller (2001) was stained with Alexa488. The protein Paxillin has recently been shown to be separated from the substrate matrix by a distance of around 20-40nm using a 3D version of a technique based on Fluorescence Photoactivatable Localization Microscopy (FPALM) -known as interference-PALM or iPALM Kanchanawong (2010). In our study similar cells were used and the separation from the substrate surface can be expected to be comparable. In figure 2 we show a confocal fluorescence image of the cell we subsequently analyze (excitation wavelength = 488nm, emission filter ≈500-550nm). We focus on a region marked by the red box. A spectral analysis of the signal in this region shows an enhancement at longer wavelengths as expected. Using the procedure outlined in the previous sections we are able to calculate zn for the pixels in this box over the range z = 1...130.T h e fluorophore concentration (rho) was initially taken to be constant up to 2.4μma n dt h e presented histogram in figure 2 is the result following four iterations of redefining rho according to the calculated zn in a n=5 parameter space (see previous section). The size of the bin width for the spectral analysis was 9.7nm. As can be seen in figure 2 where the distance distribution (zn is plotted as a function of z), the fluorophores and hence the Paxillin are indeed localized at around the expected distance from the substrate (20-40nm). We note that compared to the iPALM setup used for these studies our technique is, once optimized, very simple and straight forward to implement. The analysis is sufficiently fast to allow for fast real-time studies (which is often impossible with e.g. iPALM) and also more versatile since there is no requirement for special photoactivatable fluorophorescent markers. We are currently using it to study the dynamics of other labeled proteins found in the vicinity of cell-substrate adhesion sites including the zinc-binding phosphoprotein Zyxin and the microfilament associated protein VA S P . Fig. 2. Left: Axial distance profile of the protein paxillin (zn versus z) as inferred from ∆λ =9.7nm resolution spectrum (λ =480→800nm) of pixels at adhesion site. Right: image of cell. Red square shows pixels over which the spectrum was analyzed to obtain image on the left. [NIH 3T3 fibroblast with Alexa488 stained paxillin, scanning confocal microscope image at λ =519(±9.7nm), with 1.4 NA objective and λ =465nm excitation. Scale-bar = 10μm].

Modifications & extensions
A natural extension to the presented nano-sectioning method would be to combine it with a different technique which allows for lateral (xy) superresolution. One such technique that is relatively straight forward to implement is FPALM which was already alluded to above. In FPALM photoactivatable fluorophores are stochastically turned on by a weak activations source so that individual molecules separated by distances larger than the diffraction limit can consecutively be localized provided one has knowledge of their point spread function and the instrument and detector response. Numerous algorithms exist for such 2D localization purposes most of which rely on fitting a normal 2D (Gaussian) intensity distribution to a region of interest surrounding each molecule. Several such routines have been made available in the form of MATLAB codes. For densely labelled samples it is often desirable to be able to fit multiple overlapping Gaussians so that the isolated single-molecule activation/excitation condition can be slightly relaxed. This is particularly useful where even a very weak activation signal is sufficient to activate a significant portion of fluorophores -which we have found to be the case for photoactivatable fluorophores near the metallic structures of interest. A MATLAB code capable of fitting "overlapping intensity distributions" may be written employing a Levenberg-Marquardt based multi Gaussian least squares optimization algorithm. The result of a multiple Gaussian fit is demonstrated in figure 3 for the distribution from two nearby molecules. We note that such algorithms are also conducive to parallel processing and thus in principle fast dynamic studies. Considering the fact that typical amount of data needed to be processed for each microscopy experiment is in the order of several hundreds of megabytes; the possibility to use acceleration power of modern GPUs to enable near real-time observation of reconstruction is welcomed. At this moment MATLAB allows us to develop our implementations for several image processing and optimization routines and possibly validate the results from our experimental parallel implementations. In following text we will focus on the use of Image processing and Optimization toolboxes to outline our framework for detection regions of interest around sparsely activated fluorophores. We expect that reader is familiar with basics of image processing and numerical optimization methods and so we only propose possible solution for stated problem. For further details about mathematical background regarding image analysis we suggest Sonka et al. (2001). Current work is underway in our laboratory on combining this with the optical nanosectioning technique outlined above. A current obstacle appears to be the difficulty in rapidly photo-bleaching of the emitters -which is an essential feature of performing fast FPALM -near the substrates. This is a well known feature of fluorophores near metallic substrates and is understood to be due to the relative suppression of the triplet decay path, which is associated with the majority of the photobleaching. FPALM based nanosectioning as mentioned above is however still possible but at the price of lowered signal to noise ratio, and with properly designed image segmentation we are able extract regions aimed for further processing. The code shown below represents the necessary steps of the image processing required for detection of relevant fluorophores. This is based on dilation functions from the mathematical morphology toolset and subsequent segmentation with experimentally adjusted thresholding values.
Listing 2. Fluorophore segmentation based on thresholding and morphology 1 img = lsm_read(fileName); 2 img_dbl = im2double(img); 3 se = strel("disk" , 5); 4 img_dil = imdilate(img_dbl,se); 5 img_dil_bw = im2bw(img_dil,threshold); Depending on the sample studied, one may observe an uneven background signal which systematically introduce errors into the model fitting which may prove detrimental. For these reasons we estimate the background intensity map based on uniform filtering over large region, e.g.: Listing 3. Fluorophore segmentation with background correction 1 corrected_bg = img_dbl -imfilter(img_dbl,fspecial("average",50)); 2 img_corr_bg_closed_bw = im2bw(corrected_bg,0.60); 3 img_corr_bg_closed_bw_dil = imdilate(img_corr_bg_closed_bw,strel("disk",2)); The main algorithm framework for detection and precise position localization of activated fluorophores consists of two parts. Lines (1)-(9) set the size of the region of interest around the detected fluorophore candidates, upscale f actor determine up to a certain level precision for the fitting (with higher values at the expense of increased computational resources). The following lines are for loading of the acquired data-stack and the rescaling to the MATLAB native double format (images with scales from 0 to 1). The subsequent section is based on the code presented in listings 4 or 3 to detect local maxima and translate them to (x,y) coordinates.
Here the for loop iterates over the range of all (x,y) fluorophore coordinates, and the code is built from 2 parts, extraction of the image region of interest around the position of the fluorophore and fitting to the image data model (a Gaussian function based on a nonlinear iterative least-squares fitting). For the main algorithm of Levenberg-Marquard, interested readers are referred to Marquardt (1963)

Concluding remarks
In summary we have outlined the data analysis protocol for a novel fluorescent imaging technique capable of performing sub-wavelength scale optical sectioning in the vicinity of metamaterial coated substrates -MeSuMo fluorescence microscopy. An outline of the approach one would take for analyzing the response properties of particular materials and how to go about fitting data to the model is explained from a perspective that may be readily implemented into MATLAB routines. Whilst the calculation can be performed by most software and self written codes, the intuitive matrix approach of MATLAB makes it well suited for the required tasks. Certain toolboxes may prove useful, 3 however no toolboxes are essential or required. The presented technique is likely to find a range of applications in the lifesciences related to studying membrane traffic or indeed any events that can be made to occur near a coated substrate and nanoscale axial dynamics are of interest. In practice a fundamental limitations of the technique is the rate and accuracy with which the spectral information over a give region can be acquired. Whilst the use of multiple photodetectors combined with fast scanning of the sample may be the most accurate and best suited for not so bright samples, one can in principle perform the spectral separation optically using dichroic mirrors & filters and an array of sensitive CCD cameras (or several with split chips). The latter would remove the need for scanning, with the imaging rate now limited primarily by photon statistics.