Biobeam—Multiplexed wave-optical simulations of light-sheet microscopy

Sample-induced image-degradation remains an intricate wave-optical problem in light-sheet microscopy. Here we present biobeam, an open-source software package that enables simulation of operational light-sheet microscopes by combining data from 105–106 multiplexed and GPU-accelerated point-spread-function calculations. The wave-optical nature of these simulations leads to the faithful reproduction of spatially varying aberrations, diffraction artifacts, geometric image distortions, adaptive optics, and emergent wave-optical phenomena, and renders image-formation in light-sheet microscopy computationally tractable.


Supplementary
: Validation of biobeam with analytical solutions. a) Plane wave scattered by three solid spheres (λ=500nm, r=2-2.5µm, refractive index contrast m=1.05), b) Comparison of analytical solution (Mie calculus) versus biobeam simulation. c) Error percentage of near field distribution as a function of single sphere radius r (∆n = 0.05) and refractive index contrast ∆n (r=2.5µm). d) Top: Phase function of analytically tractable coated spheres as cell models (m=1.02/1.04, r=5µm/4µm) shows high accuracy up to approximately 0.5 radians. Bottom: size dependent scattering efficiency of the same sphere architecture and its inverse.
Supplementary Figure 2: Propagation of different predefined input fields through a tissue model of size (100µm, 100µm, 100µm) and grid dimension (1024 3 ). The respective pupil function is shown in the upper row.
Supplementary Figure 3: Detection aberration and PSF calculation. Propagating a diffraction limited input field through parts of the sample and refocusing by an idealized optical system gives the focus field as seen by the detector. If the refocus spots are separated for different starting points, the propagation of a complete grid can be carried out in a highly multiplexed manner, accelerating the process for typical microscopy simulations by a factor 100-1000.  with an incoherent light source (M470L3 Thorlabs, λ 0 = 470nm) such that an almost plane wave (NA = 0.001) illuminated the knife edge. The diffracting light was imaged below at different depths from the edge. b) The simulation was done on a computational cell of size (1024 × 256 × 1830) with voxel size ∆x = 0.29µm. We simulated the diffraction in the case of a single plane wave (coherent, top) and the incoherent superposition of 100 incident plane waves of uniformly sampled wavelengths λ ∈ [460nm, 480nm], corresponding to the measured spectral width of ±10nm of the light source (incoherent, bottom). c) The experimentally acquired intensity. Scale bar is 12µm in both axial and lateral direction (depicted with axial/lateral aspect ratio of 8, due to space constraints). d) Intensity plot at a given axial position (dashed line) for simulation, experiment and the intensity calculated via Fresnel-integral (Theory). b c simulation experiment schematic a Supplementary Figure 6: Experimental validation on a commercial light-sheet microscope. a) Polymethylmethacrylate (PMMA) microparticles with a diameter of 20µm and refractive index of n = 1.495 were embedded in an block of OptiPrep (Progen Biotechnik GmbH) / agarose (Sigma Aldrich) with refractive index of n ≈ 1.43) and which was labelled with Alexa Fluor 488. A stationary illuminating light sheet with a waist of 1.7µm and a lateral extension of ≈ 100µm was generated with a LZ1 (Zeiss) light-sheet microscope, incident on the agarose embedded sphere. Stacks were acquired at a step size of 0.414µm. b) Simulation results of the intensity distribution behind the sphere at a plane incident to the sphere center. c) Experimental intensity image. Scale bar is 20µm in both cases. Dashed lines indicate regions with specific diffraction patterns that the simulation correctly reconstitutes.  Figure 10: a) A rectilinear stripe pattern occurs distorted when imaged through a sphere. b) Fitting a radial distortion map to simulation results of this scenario, allows one to fix the aberrations at the region of interest (green) as seen through the spherical cell phantom.

Supplementary Notes
1 Numerical methods 9 1 Numerical methods

51
At the lowest level biobeam currently uses the well described scalar beam propagation method (BPM

52
[3, 4]) along with locally reduced refractive index contrasts and a mathematically exact propagator 53 whose use we explain and justify in the following:

54
The simulation of light propagation through tissue amounts to solving for the electrical field E(x, y, z) given a refractive index distribution n(x, y, z) and certain boundary conditions. This in general requires the numerical treatment of the time dependent vectorial Maxwell's equations [5]. For monochromatic illumination along z and low refractive index differences however, a far simpler description in terms of a complex scalar field u(x, y, z) becomes applicable, and the problem reduces to solving the scalar y is the accurate (i.e. non paraxial) propagator in the Fourier domain. In a further approximation the refractive index is assumed to be a small variation around a constant n(x, y, z) = n 0 (z) + ∆n(x, y, z) so that the final approximation gives the scalar beam propagation which can be efficiently solved by operator split stepping via FFTs and point wise multiplications [3,

62
To simulate the actual light propagation, we first assume the tissue to be given by a complex grid of   Crucial for the performance of BPM is a good choice of the refractive index representation n 0 (z) while 68 propagating the fields because BPM otherwise quickly becomes inaccurate for non-paraxial light fields.

69
The problems one may encounter are twofold: i) Phase shifts at large angles are under-represented 70 when refractive index contrasts are high, since the increasing path length with angle is not represented. Apart from its technical focus on speed, biobeam is specifically designed to make wave optical experi-106 ments in-silico as easy as possible. As an example, the listing 2 shows how to propagate a Bessel beam 107 with the annulus defined by the apertures NA 1 = 0.4, NA 2 = 0.43 and focal point 50µm through a 108 random refractive index volume (n = 1.33 ± 0.05) of size (100µm) 3 on a grid volume (512, 512, 512) 109 and returning the complete complex field on the grid.

157
We further calculate the scattering phase function f (θ, φ) and scattering cross section σ from the angular spectrumŨ (k x , k y , z last ) of the last plane via f (θ, φ) = −ik 0 cos θŨ (k 0 cos φ sin θ, k 0 sin φ sin θ) For a concentrically coated sphere as an analytically tractable cell model, with a cytoplasmic fraction 158 and higher refractive cell nucleus, the results are near identical (Supp. Fig.1 c) in both near and far field.
For each constant slice at z the intensityf (x, y, z) gets convolved with h(x, y, z F − z) resulting in and the overall observed image plane g(x, y, z F ) for a fixed focal position is then the integral of all Note that if the illumination field can be factored as j(x, y, z) = j xy (x, y) · j z (z), then g(x, y, z F ) = (h · j z ) ⊗ (f · j xy ) (8.5) and the observed image g is simply the 3 dimensional convolution of the illuminated density f · j xy with 268 an effective PSF h ef f = h · j z . For a distorted illumination field, this however is not true in general and significantly overlap, hence can be propagated and refocused in multiplexed way (cf. Supp. Fig. 3).

291
The PSF grid sampling is practically done at a grid-spacing below the smallest isoplanatic patch size In practice, it might often be desired not only simulate sets of PSFs distributed over a sample, but to 304 directly extract spatial maps of aberrations that would result when imaging a tissue or embryo that is 305 represented in a refractive index model.

306
Given the efficiently multiplexed determination of PSF as described in the previous section (Supp. Fig. 3), also aberrations terms can efficiently be calculated for any point in a sample in the same manner.

308
The main benefit again is that a single multiplexed simulation through the tissue is sufficient. Let 309 h(x − x 0 , y − y 0 , z F ) be the lateral slice of the PSF associated to a given initial position (x 0 , y 0 ) at focal 310 position z F , then the corresponding pupil function P (θ, φ) is given by the properly rescaled Fourier 311 transform of h(x − x 0 , y − y 0 , z F ). Projecting P (θ, φ) onto a Zernike basis yields the associated Zernike 312 terms of the aberrations present at (x 0 , y 0 , z F ). Supp. Fig. 4 shows a demonstration of these aberration 313 calculations for a tissue model with a biological plausible refractive index distribution mimicking cell 314 nuclei, cytoplasm and the eggshell of of n ∈ (1.35, 1.43) [1].

315
Biobeam is highly efficient in determining sample induced aberration for given refractive index models 316 (cf. Supp. Fig. 4). For this, light from diffraction limited guide-stars is propagated to the back-pupil 317 of a virtual microscope. Recorded field can then be composed according to Zernike aberration modes.

318
For this we record the complex scalar product between the field distribution and each Zernike mode 319 on a pupil representing unit-sphere.

329
To simulate the memory effect as in Fig. 2, we first created synthetic refractive index distributions 330 at different depth (N x /N y = 100µm, N z = 20 . . . 80µm) and with refractive index variation of n = 331 1.35 ± 0.03 by either randomly placing hard spheres in the volume or using generated Perlin noise of 332 the given variation. Next, diffraction limited guide stars were implemented inside this tissue dummy 333 and propagated to its surface, where the guide-star specific aberration patterns φ g were recorded.

334
Phase conjugation of these fields at the tissue surfaces leads to the precise recovery of the initial PSF 335 inside the tissue. Reduction of Strehl ratio was then determined as the average intensity of the focus 336 that is recreated for laterally displaced fields at the tissue surface. The decay of focal intensity with 337 lateral shifts were found identical to the absolute value of the spatial correlation in a speckle field 338 that results from an incident plane-wave. The area of the isoplanatic patch was calculated as the 339 radial integral over the Strehl ratio increase gained from the aberration pre-compensation. This is