A Computational Framework for Bioimaging Simulation

Using bioimaging technology, biologists have attempted to identify and document analytical interpretations that underlie biological phenomena in biological cells. Theoretical biology aims at distilling those interpretations into knowledge in the mathematical form of biochemical reaction networks and understanding how higher level functions emerge from the combined action of biomolecules. However, there still remain formidable challenges in bridging the gap between bioimaging and mathematical modeling. Generally, measurements using fluorescence microscopy systems are influenced by systematic effects that arise from stochastic nature of biological cells, the imaging apparatus, and optical physics. Such systematic effects are always present in all bioimaging systems and hinder quantitative comparison between the cell model and bioimages. Computational tools for such a comparison are still unavailable. Thus, in this work, we present a computational framework for handling the parameters of the cell models and the optical physics governing bioimaging systems. Simulation using this framework can generate digital images of cell simulation results after accounting for the systematic effects. We then demonstrate that such a framework enables comparison at the level of photon-counting units.


Implementation
We proposed a standard computational framework to simulate the passage of photons through fluorescent molecules and the optical system, and to generate a digital image that closely represents the image obtained using an actual fluorescence microscopy system. The computational framework included a statistical model of the systematic effects that are influenced by the parameters defined in the cell model and optical system. Using this framework, we specifically implemented the simulation modules for relatively simple bioimaging systems: TIRFM and LSCM techniques. Optical configurations are shown in Figures S1.1 and S1. 22. Those modules were designed to generate digital images that closely represent the actual digital images obtained using actual TIRFM and LSCM systems. The optics simulation of the passage of photon through the microscopy systems was based on either geometric optics or the Monte Carlo method. The optics simulation is composed of three components; (1) illumination system, (2) molecular fluorescence and (3) the image-forming system. The illumination system transfers the photon flux from a light source to a cell model, to generate a prescribed photon distribution and maximize the flux delivered to the cell model. Fluorophores defined in the cell model absorb photons from the distribution, and are quantum-mechanically excited to higher energy states. Molecular fluorescence is the result of physical and chemical processes in which the fluorophores emit photons in the excited state. Finally, the image-forming system relays a nearly exact image of the cell model to the light-sensitive detector.
Implementations of the simulation modules are generally not practical and require much time. Assuming the first-order paraxial approximation and the spatial PSF integration to be unity within a limited volume region ( ∫ Λ 0 P SF d 3 r = 1), we implement the TIRFM and LSCM simulation modules. Theoretical formulas are often applied in the first implementation. Simulation studies to estimate the errors that arise from the PSF normalization and optical aberrations are required for future implementation.

TIRFM simulation module
The TIRFM simulation module enables a selective visualization of basal surface regions of a cell model. Its optical configuration is shown in Figure S1.1 [35,33]. Implementation assumption are summarized in Table S1.1.

Fluorescence
Linear convertion (×10 −6 ) Cross-section (σ ∼ = 10 −14 cm 2 ) Image-forming 3-D Airy PSF (Unpolarized analytical form) CMOS or EMCCD cameras S1 An incident beam of excitation wavelength (λ) that passed through the objective lens is assumed to uniformly illuminate specimen. The survived photons through the use of excitation filters interact with the fluorophores in the cell model, and excite the fluorophores to the electrically excited state. The optics simulations for the focusing of the incident photons through the objective lens include a statistical model of systematic parameter rules by specifications including numerical aperture (NA), magnification, working distance, degree of aberration, correction of refracting surface radius, thickness, refractive index and details of each lens element. Details of the illumination optics are described in refs. [16,17].
The incidence angle of the beam is particularly important for the TIRFM system. Figure S1.2 illustrates schematic views of the evanescent field and epifluorescence fields with respect to different incidence beam angles. If the incidence angles are less than the critical angles given by sin θ c = n 2 /n 1 , then most of the incidence beam propagates through the interface into the lower index material with a refraction angle given by Snell's Law. However, if the incidence angle is θ > θ c , then the incidence beam undergoes total internal refraction (TIR). The evanescent field is generated along the z-axis as perpendicular to the TIR surface, and is capable of exciting the fluorescent molecules near the surface. The intensity of the evanescent field at any position exponentially decays with z, and is written in the form of where E T and A T are the transmitted electric field and amplitude of the incident beam as a function of incident beam angle, respectively. d and λ are the penetration depth of the evanescent field and the wavelength of the incident beam in vacuum, respectively. The polarization of the evanescent field depends on the incident beam polarization, which can be either p-pol (polarized in the plane of the incidence formed by the incident and reflected rays, denoted here as the x-z plane) or s-pol (polarized normal to the plane of incidence; here, in the y-direction). In both polarizations, the evanescent field fronts travel parallel to the surface in the x-direction. The p-pol evanescent field is a mixture of the transverse (z) and longitudinal (x) components; this distinguishes the p-pol evanescent field from freely propagating subcritical refracted light, which has no component longitudinal to the direction of travel. The longitudinal x component of the p-pol evanescent field range diminishes back toward the critical angle.
For the s-pol evanescent field, the evanescent electric field vector direction remains purely normal to the plane of incidence (y-direction).
where AI is the field amplitude of the polarized incident beam. The phases relative to the incident beam are written as follows.
The incident electric field amplitude in the substrate is normalized to unity for each polarization. More details are described in refs. [32,33,34,35,36].

Molecular fluorescence
Incident photons of specific wavelengths are absorbed by the fluorophores in the cell model. Fluorescence is the result of physical and chemical processes in which the fluorophores emit photons in the electronically excited state. The excitation of the fluorophores by an incident beam photon occurs in femtoseconds. Vibrational relaxation of excited-state electrons to the lowest energy state is much slower and can be measured in picoseconds. The fluorescence process (that is, emission of a longer-wavelength photon and the return of the molecule to the ground state) occurs in a relatively long time of nanoseconds. However, the process of phosphorescence from the triplet state and back to the ground state occurs in a much longer time of microseconds. A Jablonski diagram of the fluorescence process is shown in Figure S1.3. The Monte Carlo simulation of the overall fluorecence process includes a statistical model of systematic parameters ruled by the observable changes in absorption and emission spectra, quantum yield, lifetime, quenching, photobleaching and blinking, anisotropy, energy transfer, solvent effects, diffusion, complex formation, and a host of environmental variables. Details are described in refs. [18,19]. S1 Fig. S1.3: Jablonski diagram of molecular fluorescence and phosphorescence

Image-forming system
In an optical system that employs incoherent illumination of the cell model, the image-forming process can be considered as a linear system. The impulse response of an image-forming system to a pointlike fluorophore is described by the point spread function (PSF) of wavelength and position. The optics simulations of PSF formation and convolution include a statistical model of systematic parameters ruled by the observational changes in specifications of the objective lens and special filters. When all the fluorophores in the cell model are imaged simultaneously, the distribution of emitted photons of longer wavelengths passing through the objective lens and special filters used is computed as the sum of the PSFs of all the fluorophores. In the TIRFM system, the incident light that excites the fluorescent molecules is an evanescent field generated under the total internal reflection. The polarization of this light is non-isotropic, which means that dipoles of different orientations are excited with different probabilities per unit time. Therefore, the PSF of a fluorescent molecule should be written in the polarized form of the weighted average over orientations. Here, we use a simple analytical form of the unpolarized PSF models. The number of outgoing photons can spread depending on the PSF. The analytical forms of the PSF models are well studied in ref. [22] and is generally written as n . λ ′ and n are fluorophore wavelength and the refractive index of the immersion layer. The phase factor, ψ = ψ(r, z, ρ) enables generating the second Airy peak along the z-axis. In particular, we use Born-Wolf PSF model given by and C is the normalisation factor. In this PSF model, when the particle is in focus, the scalar-based diffraction can occur in the microscopy system, but the imaging plane is not required to be in focus. The model is shift invariant in all directions. A schematic view of the model condition is shown in Figure S1 The Monte Carlo simulation of a camera system includes a statistical model of a systematic source and generates digital images that closely represents the actual image obtained using an actual camera. We particularly simulate the detection process for CMOS and EMCCD cameras. Details of the camera simulations are described as follows.
(1) Uncertainty sources: Uncertainty sources of the camera systems are ruled by camera specifications and conditions shown in Table S1.2 [4,5]. First, shot noise arises from statistical fluctuations in the number of photons incident to the camera. This noise source is a fundamental property of the quantum nature of light and is always present in imaging systems. The incident photons interact with the photodiode placed on a pixel plate. Photoelectric effects can convert the incident photon signals to photoelectrons. The probability for such a conversion is the so-called quantum efficiency (QE). As both photons and electrons are quantized, the detection process is characterized by binomial distributions. Finally, readout noise is generated whereas the photoelectron signals can be linearly digitized as an image in terms of the count 16bit analog-to-digital converter (ADC). For CMOS and EMCCD cameras, the linear relationships of photoelectrons outputs with ADC outputs are shown in Figure S1.5.
n addition, the EMCCD camera has excess noise that increases the standard deviation of the output signal by √ 2 [32,37,38,39]], whereas the CMOS and CCD cameras have no excess noise (1.0). The EMCCD camera uses the multiplication process, and each stage has a small gain to multiply the number of photoelectrons. Such process is stochastic and characterized by multistage binomial distributions, which increased noise. (2) Probability density function (PDF) per pixel: The camera pixel output is the convolution of the probability distributions of each of the systematic sources. The PDF of CMOS camera pixels is given by the Poisson distribution and written in the form of where S i and E i are the random number of output electrons and expectations in the i-th pixel, respectively. The left panel of Figure S1.6 shows the PDF with respect to the number of incident photons. The PDF of EMCCD camera pixels [32,37] is written in the form of where I 1 is the modified Bessel function of the first kind of order one. α is the inverse of the EM gain. Figure S1.7 and the right panel of Figure S1.6 show the PDF with respect to the number of incident photons. S1 Fig. S1.6: Probability density function for CMOS (left) and EMCCD (right) S1 Fig. S1.7: Probability density function for EMCCD : EM gain ×100 (left) and ×300 (right) (3) Readout noise: Noise triggered by the readout electronics is typically dominated by the noise on the floating diffusion amplifier and the A/D converter. It increases with clocking speed or frame readout speed. This noise is the result of the statistical uncertainty that occurs when the amplifier attempted to reset itself to zero before the next image. The readout noise distribution for the EMCCD camera is usually Gaussian, as shown on the left panel of Figure S1.8 and Figure  S1.9. However, the readout noise distribution for the CMOS camera is uneven, because of the differences in characteristics of the amplifiers in each pixel. The distribution is shown on the right panel of Figure S1.8 and Figure S1.10.
The input photon signal (S) and optical background (I b ) falling on the photodiodes have average photon flux per pixel. The fluctuations at this rate are governed by Poisson statistics and therefore have a standard deviation that is the square root of the number of photons ( √ S + I b ). The quantum efficiency (QE) of the camera is the wavelength dependent probability that photon is converted to a photoelectron. Since the QE is predominated in the SNR equation, high QE is a fundamental attribute for obtaining high SNR. Readout noise (N r ) is a statistical expression of the variability within the electrons that convert the charge of the photoelectrons in each pixel to the number of ADC counts. EM gain (M ) is a factor of electron multiplication. It occurs in voltage dependent, stepwise manner and the total amour is a combination of the voltage applied and number of steps in EM register. The EM gain also has a statistical distribution and an associated variance accounted for the excess noise factor (F n ). The SNR and relative SNR for three cameras specifications are shown in Figures S1.11 and S1.12. S1 Fig. S1.11: No background photons, SNR (left) and relative SNR (right) for CMOS and EMCCD S1 Fig. S1.12: Five background photons, SNR (left) and relative SNR (right) for CMOS and EMCCD Figure S1.13, S1.14 and S1. 15 Figure S1. 16. We assumed that 100 TMR molecules are stationary, and randomly distributed on the surface (30 × 30 µm 2 ). Images are simulated for the optical specification and condition of the TIRFM simulation module shown in Table S1.3. Right of Figure S1.16 is the expected image averaged by 160 images captured with CMOS camera. Results are shown in Figure S1.17. S1 Figs from top row to bottom one correspond to the beam inputs of 20, 30, 40 and 50 W/cm 2 , respectively. S1 Fig. S1. 16

CMOS EMCCD ×300
EMCCD ×500 S1 Fig. S1.17: Comparison of single molecule images (100 × 100 pixels at image center). Increasing the beam flux density results in relatively smaller noise. Grayscale is the count number of 16-bit ADC.

Simple model 2 :
Using the TIRFM simulation module, we simulated imaging the basal region of the simple cell model of TMR-tagged molecules diffusing on membrane and in cytoplasm. The images are simulated for the optical system and detector specification and conditions shown in Table  S1. 4. Results are shown in Figure S1.18. We assume the simple cell that express the molecules tagged with TMR fluorescent protein. Figure S1.18A and S1.

Simple model 3 :
Using the TIRFM simulation module, we simulated imaging the basal region of the two state model of TMR-tagged molecules diffusing on membrane. The images are simulated for the optical system and detector specification and conditions shown in Table S1.5. Results are shown in Figure S1.20. We assume the simple cell that express the molecules tagged with TMR fluorescent protein. Figure S1.20 shows geometry of the two state model (20 × 20 × 4 µm 3 ). The model consists of 1 chemical species and 2 kinetic parameters. 200 TMR-tagged molecules are distributed on the cell membrane and fast diffuse with 0.2 µm 2 /sec. 300 TMR-tagged molecules are distributed on the cell membrane and slow diffuse with 0.02 µm 2 /sec. S1 Fig. S1.20: geometry of the two state model.  Table S1.5: TIRFM specifications and condition to image the two state model.

LSCM simulation module
The LSCM simulation module enables a selective visualization of focal regions of cell model. Its optical configuration is shown in Figure S1.22. Implementation assumptions are summarized in Table  S1.6.

Illumination
Gaussian beam profile Continuous / Gaussian / Unpolarized

Illumination system
Uncertainty sources of the illumination system are ruled by specification and conditions of Gaussian beam profile. We assume ideal Gaussian laser-beam intensity profile, which corresponds to the theoretical T EM 00 mode. The Gaussian beam wavefront of excitation wavelength continuously illuminate the specimen, and propagate perfectly flat with all elements moving precisely in the parallel direction. The wavefront quickly generate the 1/e 2 irradiance contour at the plane after the wavefront has propagated a distance z. The contour spreads in the form of where w 0 is the beam waist radius at the focal plane where the wavefront is assumed to be flat. z is the distance of propagation from the focal plane where the wavefront is assumed to be flat. λ is the wavelength of excitation. The intensity of the Gaussian TEM 00 beam is written in the form of where P ′ is the beam flux at the level of photon-counting unit [#photons/sec]. More details are described in ref. [17,43].

Image-forming system
The Monte Carlo simulation for the photomultipliers tube (PMT) includes a statistical model of noise source. Emitted photons of longer wavelengths are distributed as the sum of the PSFs shown in Eq. (9). As the incident beam is scanned across the cell model in horizontal and vertical axes, a digital image that closely represents the actual confocal image can be obtained pixel-by-pixel. Details of the PMT simulation are described as follows.
(1) Uncertainty sources [17,44]: Uncertainty sources of the PMT system are ruled by PMT specifications and conditions shown in Table S4. First, shot noise arises in the number of photons incident to the PMT. When the incident photons interact with the photocathode placed on the head part of a PMT, photoelectrons are emitted. These photoelectrons are multiplied by the cascade process of secondary emission through a series of dynodes and finally reach the anode connected to an output processing circuit. The methods of readout processing the output signal of a PMT can be broadly divided into the analog and digital (photon counting) modes, depending on the number of incident photons and the bandwidth of the output processing circuit. If the output pulse-to-pulse interval is narrower than each pulse width or if the signal processing circuit is not sufficiently fast, then the actual output pulses overlap and become a direct current with shot noise fluctuations. This method is in the analog mode. In contrast, if the output pulse intervals are separated from noise pulses, discrete output pulses can then be detected by the photon counting method. (2) Probability density function (PDF) [45] : An approximate PDF at the output of the PMT is written in the form of where I 1 is the modified Bessel function of the first kind of order one. The PMT is characterised by its average gain A and the number of dynode stages ν. The variance of the PMT output is 2AB where B = 1/2(A − 1)/(A 1/ν − 1). Assuming A = 10 6 and ν = 11. E is the number of photoelectrons emitted at photocathode (expectation). Approximated formulas are described in ref. [46]. The PDFs in the photon-counting mode and analog mode are shown in Figure S1.23. S1 Fig. S1.23: Probability density function for photon-counting mode (left) and analog mode (right) (3) Count rate and linearity [44]: The photon counting mode offers excellent linearity over a wide range. The lower limit is determined by the number of dark current pulses. Maximum count rate is limited by the pair-pulse time resolution where two pulses can be separated at a minimum time interval. The measured count rate is given as where N s is the input photon flux, and δt is the pair-pulses time resolution (∼ 18 nsec). Linearity is shown in Figure S1.24. S1 Fig. S1.24: Linearity for photon-counting mode (3) SNR per pixel [44] : The variance of the PMT output is given by the sum of the variances of all noise sources. The SNR and detection limits are plotted in Figure S1. 25. The SNR in the analog mode is written in the form of where I k = eQEN p and N p is the number of incident photons/sec. I d is dark current. B is bandwidth in Hz (B = 1/(2T )) and T is the observational period. G is gain factor (∼ 10 6 ). F is the excess noise (F ≈ δ/(δ − 1)) and δ is the number of dynode stages. Detection limits (SN R = 1) as a function of bandwidth is given by the following equation.
The SNR in the photon-counting mode is written in the form of where N s = QE · N p and N d is the dark count/sec. The detection limit is also given by the following equation.

Examples of outputs
Despite the absence of the scanning process, Figure S1.26 and S1.27 show images of outputs and the graphs showing the signal intensity and noise of the horizontal line at vertical center. From the top to bottom rows in each figure, 10 3 , 10 4 , 10 5 , 10 6 , and 10 7 incident photon fluxes are expected in 80 × 80 pixel squares at the image center.

Photon-counting mode
Analog mode S1 Fig. S1.26: Image output of PMT. No dark count rate Simple model : We constructed relatively simple particle model of TMR in aqueous solution as shown in Figure S1. 28. We assumed that 19, 656 TMR molecules are distributed in the solution box (30×30×6 µm 3 ), and diffuse with 100 µm 2 /sec. Images are simulated for the optical specification and condition of the LSCM simulation module shown in Table S1.8. Results are shown in Figure S1. 29