An application of bioimage simulation: cooperative binding measurement

Cooperative binding processes in living cells represent a primary and highly general mechanism for achieving the collective behaviour that emerges at the systems-level properties of macromolecules, pathways, cells and organisms. Recent progress in robotics and automated techniques for sampling and analyzing complex biological data can reduce statistical uncertainties (or errors) in measuring cooperativity. An identification of the systematic uncertainties that arise from inaccuracy in the measurement processes, however, has been often non-empirical. Here we propose a comprehensive method to quantitatively evaluate the systematic variances computed not only from a computer simulation of live-cell imaging systems, but also image processing and pattern recognition algorithms for biological images. We then demonstrate that such non-statistical variances can affect our biological interpretation of cooperativity in a binding system of interest. Within this computational scheme, the biological interpretation can be more objectively evaluated and understood under a specific configuration of the measurement processes


Introduction
For the past century, biologists have attempted to identify and document qualitative properties that underlie molecular binding processes in living cells [1][2][3][4]. Cooperativity encompasses a wide range of the binding processes such as allostery, preorganization, protein folding, conformational change of molecular structures and cellular communication [1]. For example, oxygen binding to hemoglobin can exhibit a switch-like collective response with a threshold in a concentration range of stimuli [4]. The receptor system coupling to G-proteins can also display such collective response, but is less decisive and less restricted with respect to the range of signal molecule concentration [2].
Our increasing ability to measure the details of binding processes, in the context of single molecules, multimolecular networks and whole-cells, has a great impact on our understanding of the cooperative functionalities in a variety of binding systems: in particular, automated techniques in the statistical analysis of complex biological data have been increasing the precision of the binding measurements [5][6][7][8]. However, despite advances in computational techniques, the measurement process has overlooked quantitative evaluation of the systematic uncertainties that arise from inaccuracy in experimental data acquisition and analysis [9][10][11]. This can cause erroneous identification and interpretation of cooperativity. For example, an over binding of epidermal growth factor (EGF) ligands to EGF receptors in an initial condition can alter apparent cooperative characteristics of the EGF receptor system [12][13][14].
A key challenge to evaluating the systematic uncertainties is finding meaningful and nonintuitive variances. There exist a large number of systematic sources in the measurement processes, but have no significant influence on qualitative interpretation of the measured biological properties: an equilibrium binding curve and Scatchard plot are often used to identify cooperativity [3,4]. Although experienced experimental biologists anticipate some systematic sources and make sure that most systematic uncertainties are much less than the required precision, others cannot even be recognized. There is no experimental procedure for evaluating the systematic variances.
A computational modeling approach of a whole experimental system is more relevant to estimate the systematic variances [15][16][17]. In this article, we introduce a comprehensive method to quantitatively evaluate the systematic variances computed not only from a computer simulation of live-cell imaging systems, but also image processing and pattern recognition algorithms for biological images. In particular, we construct a simulation module for the oblique illumination fluorescence microscopy system configured to observe the receptor proteins binding to ligands on an apical surface of biological cells. We then show how the non-statistical variances lead to misidentification of cooperativity in the binding system.

Computational method
Image analysis and algorithms are generally applied to actual microscopy images (or photomicrographs) in order to extract (or reconstruct) biological properties of the binding system [18]. For example, the ATLAS spot-detection algorithm can automatically capture the optimal scale corresponding to the most frequent spot size in fluorescence microcopy images [19]. However, an adequate amount of quantitative information to interpret biological meaning of the reconstructed properties has been unprovided: the mere extraction of biological properties from actual microscopy images (or experimental data) is insufficient to support our biological interpretation of the binding system.
The key insight into evaluating the systematic variances is to estimate how well ground-true model properties can be restored through the analytical procedure, influencing the biological interpretation of the reconstructed properties [9][10][11]. More precisely, quantifying the restoration efficiency and defects (or failure) in the reconstruction processes should support our understanding of the experimental outputs that arise from various combinations of model-parameter values over space and time. We cannot truly know the truth; otherwise, there would not be any uncertainty in our biological measurements. We may expect approximately what it should be, however, either from earlier experimental results or from the biological models derived from various experimental knowledge. Such approximations can function as a guide, but we must always compute the errors to establish the validity and confidence of our experimental results in a systematic way from the experimental configurations and conditions. The systematic variances generally cause the reconstructed properties to be shifted in one direction from the ground-true model properties [9][10][11]. Of particular importance is quantifying the systematic variances that are beyond the sensitivity and limitation of experimental configuration. If the groundtrue property is well-restored through the reconstruction process, then the true property remains unchanged, weakly influencing the biological interpretation. However, if the reconstruction process poorly restores the true property, then the reconstructed property can be different from the true property, thereby affecting to the biological interpretation. For example, if such systematic variances exceed the experimental limits, the variances cannot be even detected, leading to "wrong" specification and interpretation of the reconstructed biochemical kinetic rates [11].
Realistic ground-true model simulation that approximately represents whole-experimental system can be constructed for the evaluation and interpretation. The biological measurements using live-cell imaging techniques are generally composed of various natural laws and principles of biochemistry and physics. Models of each imaging processes can be simulated within the limited ranges and dimensions of model-parameter spaces. There exist various model-simulation studies: some are physics simulations for molecular fluorescence and optical apparatus [20][21][22][23], and others are systems biology simulations aimed at explaining and predicting biological phenomenon [24][25][26]. We integrate those model-simulations into the unified model-simulation approximately corresponding to the wholeexperimental system, thus helping explore the extended dimensions of model-parameter spaces that can affect potential imaging and analytical outcomes.
Bioimage simulation is one of the multidisciplinary computational approaches for handling a large range of biochemical and physical parameters that govern image-based measurement systems [27] † . A bioimage simulation module for a specific imaging technique can be constructed to generate computational photomicrographs that arise from various systematic sources: the spatiotemporal model of biological cells, photophysics and imaging apparatus. The photomicrographs are presented at the level of photon-counting units. We have implemented the simulation modules for total internal reflection fluorescence microscopy (TIRFM) and laser-scanning confocal microscopy (LSCM). The simulation modules were designed to generate the photomicrographs that closely represent the actual digital images obtained using actual fluorescence microscopy systems. In recent years, such simulation platforms have been developed for specific imaging techniques: single-molecule microscopy, super-resolution fluorescence microscopy, fluorescence correlation spectroscopy (FCS) and fluorescence recovery after photobleaching (FRAP) [28][29][30][31][32][33][34]. Within the simulation platforms, the biological interpretation can be more objectively evaluated and understood under a specific experimental configuration of interest.

Result
The computational method presented here enables us to evaluate the systematic variances that arise from inaccuracy in the cooperative binding measurement using fluorescence microscopy. In particular, we construct a bioimage simulation module for the oblique illumination fluorescence microscopy system configured to observe biochemical reactions and aggregation of the HRG ligand induced ErbB receptors on an apical region of the cell membrane [see sections A and B in the supporting information (SI)]. We then evaluate the cooperativity in two cell-models: simple binding of the single ligands to the single receptors (see SI section C.1) and dimer formation of ligand-enhanced receptors (see SI section C.2). The result for the dimer formation is as follows.
M. Hiroshima and his colleagues have used the oblique illumination technique to acquire singlemolecule images in an equilibrium region, and then identified apparent (or observed) cooperativity to be "negative" for the dimer formation [35]. In their experiment, the ErbB receptors on a cell-surface that expresses the HRG ligands tagged with tetramethylrhodamine (TMR) fluorescent molecules appear as spots (or blobs) in the single-molecule images. Spot-detection algorithm was applied to the images not only for extracting spot-properties (e.g., intensity, size), but also reconstructing the spot area-density (spots/µm 2 ) on the cell-surface. Converting the reconstructed area-density to the equilibrium binding curve and Scatchard plot, the cooperativity was finally identified in the dimer formation.
Using the bioimage simulation module, we generated the single-molecule images simulated under the model-assumption that true-cooperativity is negative (see SI section B). We then appied the analytical procedure to the computational images in order to quantify the differences between the ground-true properties and the reconstructed ones (see SI section C.2.2).

Computational photomicrography
First, we programed the simulation module for the single-molecule experimental system using fluorescence microscopy. This simulation module can generate computational single-molecule images of the three dimensional cell-model structures represented by the dimer and higher-order oligomer (e.g., trimer, tetramer, pentamer) formations on the cell-surface [35]. There are two major simulation components; the optical system and the cell-model.
(1) Optical system (see SI section A): The simulation for the optical system is composed of three parts; (a) An illumination system transfers the photon flux from a light source to the cell-model, to generate a prescribed photon distribution and maximize the flux delivered to the cell-model. An incident beam of excitation wavelength that passed through the objective lens is assumed to uniformly illuminate specimen. In particular, the optics simulation has been implemented for a selective visualization of apical and basal surface regions of the cell-model. Various illuminations can be configured: epifluorescence illumination, oblique illumination and total internal reflection fluorescence illumination. (b) Fluorophores defined in the cell-model absorb photons from the photon distribution, and are quantum-mechanically excited to higher energy states. Molecular fluorescence is the result of photophysical and biochemical processes in which the fluorophores emit photons in the excited state. In particular, Monte Carlo simulations have been implemented for the fluorescence processes of the Beer-Lambertz law and photobleaching/photoblinking effects. (c) Finally, an image-forming system relays a nearly exact image of the cell-model to the light-sensitive detector. In particular, we have programmed for the optics simulation for the formation and convolution of point spreading function (PSF), and the Monte Carlo simulation of the electron-multiplying charge-coupled device (EMCCD) camera detection process.
(2) Cell-model (see SI section B): For the dimer formation, the Spatiocyte cell simulation method can provide the spatiotemporal cell-model of biological fluctuation that arises from stochastic changes in the cell surface geometry, the number of ErbB receptors, HRG ligand binding and unbinding, receptor states (e.g., monomer and dimer), and the translational diffusion of each receptors [36,37]. In particular, negative cooperativity has been incorporated in the model by introducing an intermediate state transition. Figure 1A illustrates the biochemical reaction network of the dimer formation. Values of model parameters are shown in Table S5.
In addition, we extended the dimerization model to higher-order oligomer formation: trimer, tetramer and pentamer. In this trial extension, we introduced four additional equilibrium constants for the molecular reactions: two interactive constants are K 7 and K 8 that characterize respectively the trimer formation of monomer and dimer, and the tetramer formation of two dimers, and others are K 9 and K 10 that characterize respectively the (n + 1)-th and (n + 2)-th order of oligomer formations (3 ≤ n ≤ 12). The higher-order oligomers were assumed to be immobile on the cell membrane. The right panel of Figure S21 illustrates the biochemical reaction network of the higher-order oligomer formation. Values of additional model parameters are shown in Table S5.
Figures 1B represent optical arrangement and cell-surface geometry of the cell-model illuminated with an incident beam angle less than the critical angle. Microscopy specifications and the operating condition are shown in Table S6. Single-molecule images obtained by the simulation of oblique illumination microscopy can be accordingly compared with actual single-molecule images at the level of photon-counting unit. The image comparison is shown in Figures 2. A movie for the image comparison is available online at https://youtu.be/Q_RkZ3E4fAM. However, the image comparisons are insufficient to check the validity of the microscopy simulation module. A more elaborate set of verification and calibration is required in the future.

Analysis
The Laplacian of the Gaussian (LoG) method [38] was applied to detect spot-like features in the computational single-molecule images. We then counted the number of the detected spots for various observational area-cuts. Each area-cuts are fixed at an image center. The spot area-density is reconstructed for the area-size varying from 10 µm 2 to 1000 µm 2 . For 10 cell-samples stimulated to the 1.0 nM ligand, Figures 3A and B show the comparison of ground-true spot area-density (2.372 spots/µm 2 ) to the reconstructed one. Figure 3A shows the efficiency of the density reconstruction for various area-cuts. While the reconstruction efficiency is limited to 75% below 300 µm 2 area-cut, the efficiency is underestimated above the area-cut due to including the defocused cell-surface regions. Figure 3B shows the fractional occupancy of defected-spots for the various area-cuts. The detected spots are either one of true-molecule spots or the defects that arise from inaccuracy in the spot-detection algorithm. Such defects are one of the systematic sources that can misidentify the true-molecule spots. For the specific area-cut (28 spots/µm 2 ), Figure 3C shows the direct comparison of the reconstructed area-density varying in time to the actual one. Binding sites are saturated near 1.7 spots/µm 2 .
Properties of the detected spots are characterized as a Gaussian function of six parameters: spot pulse-height or intensity (N 0 ), spot-positions (x 0 , y 0 ), spot pulse-widths (σ x , σ y ) and background pulse-height (b 0 ). All detected spots are fitted to the Gaussian function. For the 10 cell-samples stimulated to the 1.0 nM ligand and the 28 spots/µm 2 area-cut, Figure 4A shows the difference between the ground-true spot-position (x true ), we confirmed that a peak is located near zero on each axes, and nearly formed as Gaussian function. Root mean squared (RMS) value represents the positional resolution to 0.78 pixels (58 nm). However, the tails of each distribution appear to be asymmetric. One of the most probable explanations for the asymmetry is because of the z-axis. In our analysis, we assumed that spots are characterized as a 2-dimensional Gaussian function, ignoring the z-axis. The 3-diemnsional Gaussian fitting may be able to resolve the asymmetry of each distribution.
Figures 4B to E show the direct comparison between the spot-properties reconstructed from the simulated and the actual single-molecule images . So far, no significant discrepancy is found in the the spot pulse-height (N 0 ), background pulse-height (b 0 ) and signal-to-noise (SNR) distributions. However, we found a large difference in the comparison of the spot-size (σ x ) distribution. While the averaged size of the reconstructed spots is about 1.35 pixels (89.63 nm), the actual spot-size is relatively large: 2.00 pixels (133 nm). There exist many systematic sources that cause such discrepancy: lightscattering, autofluorescence, and optical aberration. More detailed implementation of the simulation module is required to resolve the discrepancies. Furthermore, such direct comparisons for the spotproperties are insufficient to check the validity of the microscopy simulation module. A more elaborate set of verification and calibration is required in the future.
Finally, we reconstructed the equilibrium binding curve and Scatchard plot for 120 cell-samples stimulated with HRG ligands in the concentration range of 1.0 pM to 4.0 nM. Cooperativity generally appears in the shape and concavity of the binding curve and Scatchard plot. Figures 5A and B show the binding curve and the ratio of the reconstructed curve to ground-true curve. These curves clearly show that the shape of the true binding curve is partially restored in the reconstructed one. While reconstruction efficiency is about 80% and steady at relatively higher ligand concentration inputs (> 0.2 nM), the efficiency significantly reduces at lower ligand concentration region (< 0.2 nM). Figure3C shows fractional occupancy of the defected-spots. Approximately 13% of the detected spots are defected.
In a standard approach of biological science, the Hill function can be fitted to the binding curves to quantify cooperative characteristics in the ligand-receptor binding system. The Hill function can be written in the form of where L, B 0 , K A and n represent ligand concentration, maximum area-density of ligand binding, ligands occupying half of the binding sites and the Hill coefficient. The fitting results are shown in Table 1. The ground-true curve of the binding system exhibits negative cooperativity (n < 1), but positive cooperativity (n > 1) appears in the reconstructed curves.
The systematic shift clearly appears in the Scatchard plot. Figure 5D shows the comparison of the true Scatchard plot to the reconstructed one. The plot shows that the concavity of the true plot is not  Table 1: Results of fitting to the Hill equation. The best fit values and statistical errors of each parameters are listed. The ground-true curve and the reconstructed one exhibit negative (n < 1) and positive cooperativity (n > 1).
well-restored in the reconstructed one. While the model-truth is well-characterized as a concave-up curve that represents negative cooperativity, the reconstructed one exhibits a concave-down curve (or a straight line) that represents positive cooperativity (or no cooperativity).
What causes the systematic shift?
There are three major sources that can generate such systematic shifts: (1) The stochastic process of photon-detection is one of the systematic sources that may influence on the reconstruction procedures, changing the shape and concavity of the equilibrium binding curve and Scatchard plot. We generate the stochastic and expected images by turning on and off the noise channels in the simulation module. We then analyze those images to evaluate the systematic influence to the biological properties. Figures  6A and B show that all reconstructed properties are aligned in parallel, finding no significant change in the shape and concavity. (2) Imperfect performance of the spot-detection algorithm may influence on our biological interpretation of cooperativity. Figures 6C and D show examples of the time course data for the ligand concentration inputs more than 0.2 nM. In this concentration range, the reconstructed responses (presented in the red cross in the Figures) successfully converge to the 80% of the trueequilibrium state (pink dashed lines), reaching the 80% restoration of the true-equilibrium state. As we showed in Figure 5B, the shape and concavity of the reconstructed curves are unchanged from the ground-truth in the concentration range, implying weak influence on the identification of cooperativity.
(3) Quasi-static responses may be a critical issue for the identification of cooperativity. This process is a slow transition of the binding system to come to an equilibrium state. If the binding system is quite sensitive to such slow response, then the system cannot converge to the equilibrium state within the acquisition period. Figures 6E and F show examples of the time course data for the ligand concentration inputs less than 0.2 nM. In this concentration range, the reconstructed responses still remain in the nonequilibrium state within the image acquisition period: 0 to 5, 000 sec. Although the true-nonequilibrium binding state (black lines) is well-restored through the reconstruction procedure, the reconstructed binding state (red cross) is failed converging to the true-equilibrium states (presented in the pink lines), generating a gap between the reconstructed binding state and the true-equilibrium one. Such a gap cannot even be detected during the image acquisition period, thus leading to the misidentification of cooperativity.

Confidence in the fitting results
We estimated statistical uncertainties in the fitting parameters to indicate numerically our confidence in our fitting results [9,10]. Figure 8A shows the ∆χ 2 -contour plot for the results of the parameter fitting in the full concentration range: the ground-truth (pink point) is located out of the 3σ confidence contour line, implying the restoration failure of the true parameter values. In this analysis, the quasi-static response in the low-concentration range is the major systematic source that can generate the gap between the reconstructed binding state and the true-equilibrium state, thus causing the restoration failure of the ground-true Hill coefficient. Such a systematic gap cannot even be closed by increasing the number of cell-samples.
Alternatively, we performed the parameter fitting in the high concentration: 0.2 to 4.0 nM (see SI section C.2.3). The fitting result are shown in the last row of Table 1 and in Figures 7. Figure 8B shows the ∆χ 2 -contour plot for the fitting results: the ground-truth (pink point) is located within the 3σ confidence contour line. While the ground-truth of the receptor system exhibits the negative cooperativity, the reconstructed Hill coefficients largely fluctuate around a unity. The cooperative characteristics is thus less determinable in this analysis, but the contour line covers the ground-truth, displaying better results.

Conclusion and discussion
Single-molecule images have been often analyzed under the assumption that the binding system converges to the equilibrium state during the image acquisition period. Our result implies, however, the violation of this assumption. For proper image analysis, we suggest assuming that the system remains in the nonequilibrium states during the acquisition period. For further understanding of the systematic variances, analysis for parameter sensitivity and limitation is required in the future.
Furthermore, Hiroshima's experiment has been also configured to acquire the single-molecule images in a nonequilibrium regions [35]. In this configuration, detected spots are linked from one frame to another one. All linked spots are presented as space-time series of spots, and identified as an eventtrajectory. The physical properties of each event-trajectories are also analyzed: diffusion coefficient, event-length, event pulse-height and event-vertex. We use the bioimage simulation module to generate the computational single-molecule images in the nonequilibrium region. We then evaluate the performance of the event reconstruction procedures (see SI section C.2.4). The reconstructed events are either one of the true-molecule events or the defects that arise from inaccuracy in the event reconstruction procedures. Such defects are one of the systematic sources that can mimic true-molecule events. The results show that a relatively large fraction of the reconstructed events may be contaminated with the defects.
We presented the first example of estimating the systematic uncertainties and of reducing them to a level that allows us to draw a proper conclusion. Such systematic uncertainties are often present in image-based biological measurements, and hinder our qualitative interpretation of biological and physical properties in a complex system of interest. Quantitative evaluation of the systematic uncertainties will become more important when larger and more complex cell-model simulation become possible. In the near future, we believe that our simulation method can help bridges the gap between experiment and theory in biological sciences.    Red, yellow and blue colored distributions represent true-molecule spots for monomer (R), monomerlike dimer (rR or pR) and dimer (RR). Green represents the distribution for the defect-spots. All simulated spectra are normalized by the total number of the detected spots.

Supporting Information
This supporting information not only provides updates of new implementation for the fluorescence microscopy simulation module, but also details for the model construction and the reconstruction of physical and biological properties. A Implementation updates Figure S1: Optical configurations of the TIRFM simulation module [27].

Contents
We have implemented the simulation modules of total internal reflection fluorescence microscopy (TIRFM) and laserscanning confocal microscopy (LSCM) [27]. Those simulation modules were designed to generate digital images that closely represent the actual digital images obtained using actual fluorescence microscopy systems. In this section, we provide updates of new developments to extend the TIRFM simulation module characteristics since the first publication. There are two new features in this update. (1) Monte Carlo method is applied to simulate Beer-Lambert law and the stochastic processes of photobleaching and photoblinking effects. (2) Fluorescence microscopy enables selective illumination of the apical regions of cell-model. Epiillumination and oblique illumination techniques are implemented.

A.1 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 [27]. Implementation assumptions are summarized in Table S1. 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.

Fluorescence
Beer-Lamberts law Photophysics (photobleaching ... etc) Image-forming 3-D PSF Models (Unpolarized analytical form) EMCCD camera Table S1: Implementation assumptions for the TIRFM simulation module. The detection process for the cameras is performed with Monte Carlo simulation, where EMCCD stand for electron-multiplying charge-coupled device

A.1.1 Molecular fluorescence
Molecular fluorescence is the result of physical and chemical processes in which fluorophores emit photons from electronically excited states [20,21,23]. The Monte Carlo simulation of the overall fluorescence process includes a statistical model of the systematic effects that are influenced by the absorption and emission spectra, quantum yield, lifetime, quenching, photobleaching and blinking, anisotropy, energy transfer, solvent effect, diffusion, complex formation, and a host of environmental variables. In this study, we particularly implement for the fluorescence processes of Beer-Lambertz law, photobleaching and photobliking. Implementation details are described as follows; (1) Beer-Lambert law [20,21,23] : A fundamental aspect of molecular fluorescence is the attenuation of a photon to the properties of materials through which the photon is traveling. The left panel of Figure S3 shows the schematic view of the molecular fluorescence. The relationship of the number of photons entering volume (n 0 ) to the number of photons leaving the volume (n 1 ) is written in the form of where Φ QY , |A T | 2 , and δT are quantum yield, the transmitted beam flux density, and detection time. The absorption cross-section is given by σ = ln(10) where N A is Avogadro's number. The detector is located in a specific direction. We expect to observe the number of photons devided by an unit surface area of a sphere (4π). The amplitude of the transmitted beam flux density depends on the index of refraction, and the incident beam angle, amplitude and polarization. The absorption coefficient (A) is given by where , c and l are molar absorption coefficient (or cross-section), volume concentration and path length. In our simulation, we assume that the volume concentration, penetration depth and detection time are given by Spatiocyte voxel-volume, voxel-diameter, and time interval. Finally the expected number of photons emitted from a single fluorophore is given by In previous implementations of bioimaging simulation [27], we have assumed that the fluorescence molecules subsequently emit a single photon of longer wavelength while they absorb one million photons of excitation wavelength, and the cross-section of photon-molecule interaction is roughly 10 −14 cm 2 . Such approximation is given by The right panel of Figure   (2) Aging Mechanism [21,23,[39][40][41][42] : Statistical aging process, such as photobleaching effect, involves photodynamic interactions between excited fluorophores and molecular oxygen (O 2 ) in its triplet state, dissolved in the sample media. The aging process is governed by Levy statistics. In a simple model, we assume a limited number of photons that each fluorophore can emit (photon budget), and the presence of thermally activated barriers in the state transition between a bright state and a permanently photobleached dark state. Such transition leads to the number of photons decreasing along the time, preventing long exposure time experiments. We assume that the probability distribution associated to time scales of the photobleaching effects, is given by a power-law, and written in the form of where α and τ 0 are dimenssionless constants (0 < α < 1) and a characteristics time. In addition, the number of photons emitted before photobleaching depends on photochemical interactions of intracellular singlet oxygen cencentration and on the distance between the fluorophores and intracellular components such as protein and lipids.
on/of f are drawn at random from the probability distribution given by the equation (6).

A.1.2 Examples of single-molecule images
We constructed realistic particle model of TMR-tagged proteins on glass surface as shown in Figure  S6. We assumed that 10, 000 fluorescent molecules are stationary, and randomly distributed on the surface (30 × 30 µm 2 ). An aggregation process where many fluorophores form clusters in a colloidal suspension, is included in the model. Images are simulated for the optical specification and condition of the TIRFM simulation module shown in Table S2. Incident beam angle is 65.7 • (Evanescent fields).
Photophysical parameterizations of fluorophores are also shown in the        Table S2: Microscopy specifications and operating conditions for imaging the particle models.

A.2 Fluorescence microscopy simulation module
We implemented the simulation module for fluorescence microscopy. The simulation module enables a selective visualization of apical and basal surface regions of a cell-model. Its optical configuration is shown in Figure S1. Implementation assumptions are summarized in Table S3.

Fluorescence
Beer-Lamberts law Photophysics (photobleaching ... etc) Image-forming 3-D PSF Models (Unpolarized analytical form) EMCCD camera Table S3: Implementation assumptions for the epifluorescence microscopy. The detection process for the cameras is performed with Monte Carlo simulation.

A.2.1 Epi-illumination / oblique illumination
An incident beam of excitation wavelength (λ) that passed through the objective lens is assumed to uniformly illuminate the specimen. The surviving 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 ruled 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. Figure S9 illustrates schematic views of epi-illumination, oblique illumination and Evanescent field. The left panel of Figure S10 shows the relative intensity of incident electric fields as a function of incident beam angle for s-and p-polarization. The polarization of the illumination 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). If the incidence beam angles are less than the critical angles given by sin θ c = n 2 /n 1 where medias are fused silica (n 1 = 1.46) and water (n 2 = 1.33), then most of the incidence beam propagates through the interface into the lower index material with a refraction angle given by Snell's Law. If the incidence angle is at 0 degree, then the epi-illumination is generated along the z-axis. If the angle is at 0 < θ θ c , the oblique illumination is generated along the refracted beam axis. Finally, if the angle is θ > θ c , then the incidence beam undergoes total internal refraction (TIR). The evanescent field is generated along the z-axis, 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. The right panel of Figure S10 shows the penetration depth of the evanescent field as a function of the incident beam angle. More details are described in reference [47][48][49][50]. Examples of single-molecule images : We constructed the simple cell-model of TMR-tagged proteins distributed on membrane surface. We assumed that 10, 000 fluorescent molecules are stationary, and randomly distributed on the cell surface (30 × 30 × 10 µm 3 ). The aggregation process where many fluorophores form clusters in a colloidal suspension, is included in the model. Single-molecule images are simulated for the optical specifications and operating conditions of the fluorescence miroscopy simulation module shown in Table S2. Figure S11 shows the epi-illumination (θ = 0 • ), the oblique illumination (θ = 8 • , 16 • , 24 • ) and evanescent field (θ = 65.7 • ). Photophysical parameterizations of fluorophores are also shown in the Table. Results are shown in Figures S12-S13. Each figures show snapshot images captured with fluorescence miroscopy imaging configuration.   71st Frame (7.0 sec) Figure S13: Examples of single-molecule images (2). 8 frames of the cell-model that includes particle aggregation and photobleaching effects. Each image size is 512 × 512 pixels. Actual minimum and maximum values of the intensity histogram are 1, 900 and 4, 800 ADC counts. The intensity is rescaled in the range of 0 to 255.

A.2.2 Oblique illumination near the critical angle
Oblique illumination is designed to drastically reduce the background of scattering light originating from solution, dust and optical elements [35,[51][52][53][54][55]. This illumination technique allows us to investigate the dynamics of individual molecules at an apical cell surface. The left panel of Figure S14 illustrates schematic views of the oblique illumination near the critical angle, generating evanescent field at the apical cell surface. An incident beam of excitation wavelength (λ) that passed through the objective lens is assumed to uniformly illuminate the specimen. The surviving photons through the use of excitation filters interact with the fluorophores in the cell-model, and excite the fluorophores to the electrically excited state. If the incidence beam angles are less than the critical angles given by sin θ c = n i /n j where media are fused silica (n 1 = 1.460), cell (n 2 = 1.384) and culture medium (n 3 = 1.337) [56], then most of the incidence beam propagates through the interface into the lower index material with a refraction angle given by Snell's Law. If the incident beam angle is near the critical angle (0 θ ≤ θ c ), then the incidence beam can be reflected between cell cytoplasm and culture medium, and undergoes total internal refraction (TIR) at the apical cell surface. The evanescent field is generated along the axis as perpendicular to the apical surface, and is capable of exciting the fluorescent molecules near the apical surface. The intensity of the evanescent field at any position depends on cell geometry, and exponentially decays with the surface normal axis.
Examples of single-molecule images : We assumed a hemispherical cell-model that represents the fluorescence molecules tagged with tetramethylrhodamine (TMR). The cell measuring 50 µm in diameter, and 6 µm in height was placed on the glass surface. Approximately 8, 000 molecules are uniformly distributed on the cell surface. We then simulated single-molecule images for the optical specification and condition of the fluorescence microscopy simulation module shown in Table S4. Photophysical parameterizations of fluorophores are also shown in the Table. Photobleaching and blinking effects are not included in the example images. The right panel of Figure S14 shows relative intensity of incident electric fields as a function of incident beam angle. Results are shown in Figures S16 and S17. Each Figure shows snapshot images captured with fluorescence microscopy imaging configuration. In addition, the use of two excitation beams allows us to observe the dynamics of single-molecules for a given focal height [52]. Figure S15 llustrates schematic views of the oblique illumination using two incident beams. Assuming the simple cell-model introduced above, we simulated single-molecule images for the microscopy specifications and operating conditions shown in Table S4. Results are shown in Figures S18 and S19. Each       Figure S19: Histograms of single-molecule images (2). Two excitation beams are used to generate these images. Each image size is 512 × 512 pixels. θ 1 , θ 2 and z 0 represent the incident beam angles and the focal distance from the glass surface. Actual minimum and maximum values of the intensity histogram are 1, 900 and 4, 200 ADC counts. The intensity is rescaled in the range of 0 to 255.

B Model construction B.1 Hiroshima's model of dimer formation
We constructed a mathematical cell-model of HRG ligand induced ErbB receptor dimerization in a biological cell that represents the HRG ligands tagged with tetramethylrhodamine (TMR). A cell simulation method with Spatiocyte was used to construct the cell-model of the dimer formation. Spatiocyte provides the spatial cell-model of biological fluctuation that arises from stochastic changes in the cell surface geometry, number of ErbB receptors, HRG ligand binding, molecular states (monomer or dimer), and translational diffusion of each receptors [36]. The left panel of Figure S21 illustrates biochemical reaction network of the dimer formation. An intermediate state transition of the proteins is included in this model. The hemi-elliptical cell measuring 65 µm and 45 µm in minor and major axes, and 2.0 µm in height was placed on the glass surface. Average observational surface area of the cell-model is 4, 731 µm 2 . Monomer and dimer of ErbB receptors are presented as solid spheres where voxel radius is 20 nm. Each receptors is uniformly distributed in the surface compartment, and randomly walk from voxel to voxel, slowly diffusing at 0.015 µm 2 /sec. Molecular collisions occur between walks. We defined seven equilibrium constants for the molecular collisions: K 1 , K 2 , and K 3 are the binding constants of HRG to an isolated receptor, to a receptor in a dimer where one site is free and one bound. K 4 , K 5 , and K 6 are the dimerization constants of two free receptors, one bound and one free receptor, and two bound receptors. K i that characterizes the intermediate state transition constant in the dimer where one site is free and one bound. Values of each model parameters are shown in Table S5. We simulated the dimerization model without the systematic effects that arise from microscopy specification and its operating condition. Simulation results are shown in Figure  S20. The left and the right panels show the equilibrium binding curve and the Scatchard plot. The Scatchard plot shows a concave-up curve that represents negative cooperativity of the binding system. We consequently reproduced the properties obtained by Hiroshima's analysis [35].

B.2 Toy model of higher-order oligomer formation
In order to extend the dimerization model to higher-order oligomer formation, we defined four additional equilibrium constants for the molecular collisions: two interactive constants are K 7 and K 8 that characterize respectively the trimer formation of monomer and dimer, and the tetramer formation of two dimers; and others are K 9 and K 10 that characterize respectively the (n + 1)-th and (n + 2)-th order of oligomer formations (3 ≤ n ≤ 12). We assumed that the monomers and dimers are the subunits which polymerize on the cell membrane, and each higher-order oligomers are presented as immobile polymerized chains. Figure S21B illustrates an additional reaction network of the higher-order oligomer formation. Values of each parameter are shown in Table S5. In fact, average observational area on the cell surface (43 2 = 1, 930 µm 2 ) is larger than actual image size (34 2 µm 2 ) [35]. The observational area needed to be adjusted to 661 µm 2 . Such parameter adjustment consequently leads to changes in the total area density (1.70 → 4.96 µm −2 ) and the equilibrium constants (K 4,5,6 ).  Table S5: Model parameters for oligomer formation. Hiroshima has provided two sets of parameter values; primary and secondary parameters for the interactions between HRG and ErbB molecules [35]. We use the primary set of reaction parameter values. Monomers and dimers are properly selected for analyzing each reaction rate. However, the monomer and dimers are not properly selected for obtaining the secondary set of parameter values. The secondary set is obtained by fitting of the dimerization model to all observed molecules, including higher-order oligomers (∼ 10%).
Image comparison : We simulated single-molecule imaging of the apical surface region of each model for the optical specification and condition of the fluorescence microscopy simulation module shown in Table S6. Comparisons of actual single-molecule images to the simulated single-molecule images for 1.0 nM ligand input are shown in Figure S22-S24. The left column shows the actual images obtained by the single-molecule experiment (Movie S1 in Ref. [35]). Cell size in the actual image is relatively large, but not known. Middle and right columns show the single-molecule images obtained by the bioimaging simulations for the dimerization model and higher-order oligomerization. The photobleaching effects are included in each simulated images. The simulated images of the higher-order oligomer formation model are visually similar to the corresponding real images at steady state. Thus, the simulated images were properly compared with the images obtained using the actual fluorescence microscopy system at the level of photon-counting units. However, there still remains differences in the resulting images. A more elaborate set for calibration is required for further understanding of the receptor system.

B.3 PSF-verification
Verification is the process of confirming the simulation modules are correctly implemented with respect to conceptual description and analytical solutions [27]. During the verification process, the simulation modules must be tested to find and estimate numerical errors in the implementations. The simulation module presented here is designed to count the number of photons that passed through the optical configuration. A wrong estimation of the numerical errors that arise from the photon-counting principle can provide a wrong intensity of the final images.
In our optics simulation, an image-forming system enables the formation and convolution of the Born-Wolf form of point spreading functions (PSF). We assume that single fluorophore (voxel-radius = 20 nm, path-length = 40 nm) is placed on glass surface. We then evaluate how well the exact-from of the Born-Wolf PSF are correctly implemented for various lateral distance (z-axis) from the focal plane. In particular, two constraints are applied to optimize the PSF simulation processes: r-cut (r < 1 nm) and z-cut (z < 1 nm). Microscopy configuration is shown in Table S6. The results of the PSF verification are presented below.
We computed the total number of the emitted photons from the exact-form and the simulated-form of the PSF for two illumination configurations: oblique illumination (θ inc = 60 • ) and evanescent field (θ inc = 72 • ). The top panels of Figures  the exact-form (black lines) of the PSF with respect to the lateral distance from the focal plane. While the total number of the emitted photons from the exact-form of the PSF is constant in the oblique illumination configuration, an exponential decrease is shown in the evanescent (or TIR) configuration. The bottom panels show the ratio of the total number from the exact-form to the simulated-form of the PSF, implying that the total number is well-conserved within 4-6% differences near the focal plane (below the 600 nm lateral depth). The differences, however, increase up to 30% at defocused region or background (above the 600 nm depth). If the z-cut is omitted from the simulation processes, then the total photon number can be violated with respect to the lateral distance (red dashed lines). The additional z-cut can, however, help recovering the total photon number up to 70% levels.
We also evaluated the peak (or -maximum) photon number from the exact-from and the simulatedform of the PSF. The top panels of Figures S26 show the comparison of the PSF-peak photon number as a function of the lateral distance from the focal plane. The bottom panels show the ratio of the peak photon number from the exact-form to the simulated-form of the PSF, implying that the peak photon number is well-conserved within 2% differences below the z-cut limit. Although the difference is overflow above the limit, the peak photon number is relatively low compared with the value on the focal plane. This may weakly influence on the image-forming processes. In addition, Figures S27 and  S28 show the PSF-image comparison for various lateral distances. We found no significant difference near the focal region. The differences, however, clearly appear above the z-cut limit.

C Reconstruction of physical and biological properties
A primary task for quantitative image analysis is to accurately reconstruct physical properties of single-molecules such as molecular localization and diffusion coefficient. Our collaborators strongly rely on the reconstruction method (or single particle tracking) constructed by T. M. Watanabe in 2005 (unpublished work). In this study, we modified spot-scanning process in his reconstruction method, and computationally evaluate the performance of the modified method. There are four major steps in the modified reconstruction procedure. The steps are shortly described as follows.
(1) Spot-detection : We applied the Laplacian of the Gaussian (LoG) method to identify spots (or blobs) in single-molecule images [38] * . The method detects spot-like features by searching for scale space extrema of a scale-normalized LoG : where L(x, y; σ) = g(x, y; σ) * f (x, y) is a scale space representation at a certain scale σ. g(x, y; σ) and f (x, y) are Gaussian kernel and input image. There are five input parameters; (a) minimum and maximum standard deviation of the Gaussian kernel, (b) the number of intermediate values in the range of minimum to maximum standard deviations, (c) local maxima smaller than threshold are ignored, reducing detection of blobs with less intensities, and (d) an overlap value between 0 and 1, which can eliminate the smaller blobs if the area of two blobs overlaps by a fraction greater than the threshold. In addition, more accurate methods are discussed in ref. [19,57,58].
(2) Spot-property : Spot-properties are characterized as a Gaussian function of six parameters. Each detected spots can be fitted to the function written in the form of where N 0 , (x 0 , y 0 ), σ x,y and bg are spot pulse-height (or normalization factor), central spot position, spot pulse-widths, and background pulse-height.
(3) Spot-tracking (or event-identification) : We assume two conditions in linking two spots. The primary condition is the distance between i-th spot in k-th image-frame and j-th spot in (k + 1)th image frame must be less than average of each spot sizes. The condition can be written in the form of where i and j represent spot indexes at k-th and (k + 1)-th image-frames. The secondary condition is the intensity difference between i-th spot in k-th image-frame and j-th spot in (k + 1)-th image frame must be less than sum of each spot shot-noise. The condition is given as follows.
|n( r j,k+1 ) − n( r i,k )| < n( r i,k ) + n( r j,k+1 ) where n represents the spot pulse-height in photoelectron unit. If there is more than one spot that satisfies the two conditions, then the spot in the shortest distance can be linked. Finally, all linked spots are presented as space-time series of spots, and identified as an event.
(4) Event-property : The physical properties of each spot-series are reconstructed: diffusion coefficient, event-length, event pulse-height and event-vertex.

C.1.1 No binding
The use of simulated single-molecule images allows us to quantitatively check the performance of the spot and event reconstructions. The results are shown as follows.
Cell-model : We constructed a relatively simple cell-model where monomer and dimers slowly diffuse on a cellular membrane. The GEOMETRY and LENGTH variables of the Spatiocyte cell compartment are set to the ellipsoidal shape and 45 × 35 × 2.0 µm 3 (Y × Z × X). The radius of the HCP lattice voxels is set to 20 nm. The cell-surface area is 2576.72 µm 2 with 715 monomers and 357 dimers uniformly distributed on the cell-surface. Those molecules slowly diffuse with 0.015 µm 2 /sec. No molecular reaction is defined in this simple model.

Microscopy configuration :
Single-molecule imaging of apical and basal regions of the simple cell-model is simulated for the optical specification and its operating condition of the fluorescence microscopy simulation module shown in Table S6. Photobleaching is included in the simulated singlemolecule images. Incident beam angle is set to 72 • to observe the basal region of the cell-model. In order to observe the apical cell-area, we also set the incident beam angle to 60 • to generate oblique illumination. While exposure time is configured to 0.150 sec, image acquisition period is 10 sec (66 frames). The illumination configuration on the cell-surfaces is shown in Figure S29.  Figure S30 shows example images of spot-detection on the apical and basal cell-regions. The observational area cut to 300 × 300 pixels (397 µm 2 ) is shown with dashed lines. Red circles represent molecular-spots detected by using the LoG algorithm. Each image size is 512 × 512 pixels. Actual minimum and maximum values of the image intensity are 1, 900 and 2, 500 ADC counts. The image intensity is rescaled in the range of 0 to 255. The molecular-spots are clear and well detected in the basal cell-region. In the apical cell-region, the microscopy is configured to focused on the image center. While molecular-spots are clear and properly picked up at the focused area, the spots are blur and not well detected near the curved-region on the cellular membrane.
In addition, using the simulated images not including photobleaching effects, we reconstruct spot area-density and estimate reconstruction efficiencies for various observational cuts. The area is fixed at the image center, and the area-size varies from 10 µm 2 to 1000 µm 2 . Figure S31 shows the reconstruction efficiency of the area-density in apical and basal cell-regions. Ground-true spot areadensity is set to 0.416 spots/µm 2 . While the reconstruction efficiency in the basal cell-region is constant at 90% for the various area-cuts, the efficiency is largely varied by the area-cuts in the apical cell-regions. The reconstructed area-density is over-and under-estimated below and above 100 µm 2 area-cut.

Apical cell-region
Basal cell-region   Figure S32. This is also known as cluster size in Hiroshima's work [35]. The black solid line represents observed distribution of reconstructed spots (monomer and dimer together). Two peaks are observed in these histograms. First and second peaks correspond to monomer-like and dimer-like spot distributions. Red and blue filled histograms represent ground-truth distribution of monomer (R) and dimer (RR) spots. Some of dimer spots photobleaches and forms monomerlike spots. The histograms show that the photobleached-dimer spots (blue) are identified as  Figure S32. In those Figures, we confirmed that peaks of each resolution distribution are located near zeros, and are nearly formed as Gaussian distributions. Root mean squared (RMS) value represents the positional resolution to 0.78 pixels (58 nm). However, the tails of each distribution appear to be asymmetric. One of the possible explanations of the asymmetry is because of the z-axis. In our analysis, we assumed that spots are characterized as 2-D Gaussian distributions, ignoring the z-axis. 3-D Gaussian fitting may be able to resolve the asymmetry of each distribution.

Event-trajectories in observation
Event-trajectories in ground-truth Apical cell-region Basal cell-region Figure S35: Event-trajectories 1. Image pixel-coordinates in x-y axes is shown in the horizontal plane. The vertical axis represents time.

Trajectories for each event-types :
Examples : Apical cell-region Basal cell-region     (4) Event-property : Reconstructed track properties are presented as physical parameters: eventlength, event-vertex, diffusion-coefficient, and event pulse-height. Distributions of each physical parameters are shown in Figures S41, S42 and S43. Red, blue and green filled histograms represent the ground-truth distribution of true-monomer, true-dimer and false event-reconstruction. Black solid lines represents the the observed distributions (distribution sum of each event-types).
(a) Diffusion : The diffusion coefficient in lateral axes is reconstructed using short-range analysis method [47]. The diffusion coefficient is determined by linear-fitting to mean squared displacement (MSD) for various time intervals (δt). The fitting function is written in the form of where D and a represent diffusion coefficient and intersection in the MSD axis. Distributions of reconstructed diffusion constant and its resolution (D reco − D true ) are shown in the top panels of Figure S41. In those Figures, we confirmed that peaks of the resolution distributions are located near zeros, and are roughly formed as Gaussian distributions in logarithmic scale. The RMS value represents the diffusional resolution to 0.32 log 10 (µm 2 /sec). In addition, more cell-model samples are required to check the shape of the resolution distributions.
(b) Event-length : Time-frames and path-length are defined as the total time and distance that a molecule travels. Distributions of the time-frames and path-length are shown in the bottom panels of Figure S41.
(c) Event pulse-height : Distributions of event pulse-heights are shown in the top panels of Figure  S42. The event pulse-height is defined as the total spot pulse-height in a given trajectory, and written in the form of where N 0 is the spot pulse-height of j-th spot in a given trajectory.

C.1.2 Simple ligand-receptor binding
We construct a relatively simple cell-model of monomeric receptor binding to ligand on cellular membrane: L + r ↔ R. The monomers slowly diffuse at 0.015 µm 2 /sec. Monomer density and the dissociation constant are configured to 4.977 µm −2 and 3.63 nM (k 1 = 0.00193 nM −1 sec −1 , d 1 = 0.00700 sec −1 ). Cell-compartments are shown in Table S13. Experimental configuration is shown in Table S8.  Single-molecule imaging of the apical region of the ligand-receptor binding is simulated for the optical specification and its operating condition of the fluorescence microscopy simulation module shown in Table S6. In particular, photobleaching is not included in the simulated single-molecule images. The results are as follows.
(1) Spot-detection : In this experimental configuration, the parameter values for the spot-detection algorithm are configured to σ min = 2.0 pixels, σ max = 4.0 pixels, 20 intermediate values in the deviation range, threshold = 15, and overlap = 0.5. Figures S44 show example images of spotdetection on the apical cell-regions for 100 nM ligand concentration. We assume two observational area-cuts to 80 × 80 pixels (cut-1) and 200 × 200 pixels (cut-2) with respect to the image center.
Analysis results of various observational area-cuts are shown in Table S9 and Figures S45; (a) the Table shows the number and fraction of simulated spots. The simulated spots are true-molecular spots and false-spots that arise from molecules and background noise. In particular, the false-spots are noise-like spots that can mimic the molecular spots. Relatively large fraction of the false-spots are captured with the cut-2. (b) Top left and right Figure panels show efficiency of area-density reconstruction and fractional occupancy of false-spots for various observational cuts. The area-cut is fixed at image center, and area-size varies from 10 µm 2 to 1000 µm 2 . Ground-true spot area-density is set to 1.076 spots/µm 2 . The reconstruction efficiency is constant to 100% below 100 µm 2 area-cut. The efficiency is underestimated above the area-cut, due to including the defocused regions. (c) The left and the right bottom panels show time-course data for the cut-1 and cut-2. The reconstructed area-density varies in time, and saturated at ∼ 1.00 spots/µm 2 for each area-cuts.
(2) Spot-property : Reconstructed spot-properties are presented as the Gaussian function of six parameters; the spot pulse-height (or normalization factor), central position, spot-size (or Gaussianwidth) and background pulse-height. Distributions and correlations for each parameter is shown in Figures S46 to S48

Experimental configuration B (various ligand inputs)
In this experimental configuration, 10 cell-model samples are prepared for 12 ligand inputs (0.3 nM ∼ 400 nM). Cell-compartments for each cell-model are listed in Table S13. We also assume two observational cut-1 and cut-2. Parameter values for the spot-detection algorithm are configured to σ min = 2.0 pixels, σ max  1) The left and the right top panels show the equilibrium binding curve and its ratio of the reconstructed density to ground-truth. These plots clearly show that shape of the true binding curve is partially recovered in the reconstructed one. While the reconstruction efficiency is relatively high at low concentration, the efficiency is limited by up to 50% at relatively higher ligand concentration input. (2) The bottom left panel shows the fraction of false-spots contaminated in all detected spots. Approximately 5% of the detected spots are false-spots with the cut-1. Roughly, 65% are contaminated by the cut-2. (3) The bottom right panel shows a Scatchard plot. The plot clearly shows that shape of true-Scatchard plot is not well-recovered in the reconstructed ones. The true-Scatchard plot is well-characterized as a straight line, representing no cooperativity. The shape of the reconstructed Scatchard plot are also represented with the straight line: No significant change was found in the simple binding system.
Hill coefficient : In a standard approach of biochemistry and biology, the Hill equation (1) can be fitted to the equilibrium binding curves to evaluate signal responses in simple ligand-receptor binding. The minimization function is given in the form of where E i and O i are the expected and the observed data at the i-th bin. σ i is statistical error at the i-th bin. The fitting results are shown in Table S10. The top panels show the equilibrium binding curves for each cut. Bottom ones show the ratio of the reconstructed binding curves to the fitted ones. The Hill coefficients are nearly an unity: No significant change was found between the ground-true curve and reconstructed one, exhibiting no coopertivity in the simple binding system.  Table S10: Results of fitting to the Hill equation. The best fit values and uncertainties of each parameters are listed.χ 2 0 is the reduced minimum. If n < 1, then the receptor system increases binding affinity of states and exhibits negative cooperativity. If n > 1, then cooperativity is positive, decreasing the binding affinity of states. If n = 1, then there is no cooperativity.
Stochasticity : To see influence of stochastic photon-detection processes to the analytical procedure, we compared the biological properties reconstructed from the expected and stochastic images (see the SI section C.2.2 for more detailed explanation). Figure S51 shows the equilibrium binding curves and Scatchard plot for each area-cuts. Red and blue lines represent the reconstructed curves obtained by the cut-1 and cut-2. The reconstructed curves obtained from the stochastic images are presented with green lines. Ground-true Scatchard plots are presented with black solid lines.
The top panels show the comparison of the equilibrium binding curves reconstructed from the expected and stochastic images. The restoration efficiency and defects of the analytical procedure are shown in the middle four panels. We computed the ratio of the reconstructed binding curves of the stochastic images to the expected ones. The 2nd top panels show the ratio curves, implying that the stochasticity can increase to approximately 20% of the restoration efficiency of capturing more molecular-spots at the low concentration range (< 6 nM). The 2nd bottom panels show the fraction of false-spots contaminated in all detected spots. The left panel for the cut-1 represents that the stochasticity can increase to approximately 20% of contamination of capturing the defected spots at the low concentration range. The right panel for the cut-2, however, shows no influence of the stochasticity to the false-spots contamination. Finally, the bottom panels show the comparison of the Scatchard plots reconstructed from the expected and stochastic images. Each reconstructed Scatchard plot is crossing over, do not align in parallel. Thus the stochasticity cannot conserve the shape and concavity of the binding curve and Scatchard plot. Confidence interval : We estimate the uncertainties in the parameters to indicate numerically our confidence in our fitting results. ∆χ 2 -plots for each area-cuts are shown in Figures S52 show. Black solid lines represent the 1σ (68%), 2σ (95%) and 3σ (99%) confidence intervals. Black and pink points indicate the best fit and ground-truth. The top panels show the ∆χ 2 contour plots: the Hill coefficient (n) vs dissociation constant (K A ). The ground-truth (pink point) is located out of the 3σ confidence contour line, implying the restoration failure of the true parameter values. The middle four plots show the ∆χ 2 plots for the Hill coefficient (n) and dissociation constant (K A ).While the ground-true Hill coefficient is found within 3σ contour line, the true dissociation constant is clearly located out of the contour lines. Cooperativity is thus well-restored through the reconstruction procedure.

C.2 Dimer formation
The reconstruction procedure can be also applied for Hiroshima's model of dimer formation. Model parameters and cell-compartments are shown in Table S12 and S13.
In order to quantitatively evaluate the performance of the reconstruction procedure, two single-molecule experiments are configured at (A) the equilibrium region and (B) the nonequilibrium region [35]. Details of each experimental configurations are shown in Figure S53 and Table S11.  Table S11: Experimental configurations [35].   Table S13: Cell-compartmental variables for 10 cell-samples. GEOMETRY of each cell-model is configured to "Ellipsoid".

C.2.1 Experimental configuration A (equilibrium region)
Single-molecule imaging of apical regions of each model is simulated for the optical specification and operating condition of the fluorescence microscopy simulation module shown in Table S6. In particular, photobleaching is not included in this experimental configuration. The results are shown as follows. Analysis results of various observational area-cuts are shown in Table S14 and Figures S55; (a) the Table shows the number and fraction of observed and simulated spots. The simulated spots are true-molecular spots and false-spots that arise from molecules and background noise. In particular, the false-spots are noise-like spots that can mimic the molecular spots. A relatively large fraction of the false-spots are captured with the cut-2. (b) Top left and right Figure panels show efficiency of area-density reconstruction and fractional occupancy of false-spots for various observational cuts. The area-cut is fixed at the image center, and area-size varies from 10 µm 2 to 1000 µm 2 . Ground-true spot area-density is set to 2.372 spots/µm 2 . While the reconstruction efficiency is constant at 75% below 300 µm 2 area-cut, the efficiency is underestimated above the area-cuts including the defocused regions. (c) The left and the right bottom panels show time-course data for the cut-1 and cut-2. The reconstructed area-density varies in time, and saturated at ∼ 1.7 spots/µm 2 for each area-cuts.   Positional resolutions (or localization error) of reconstructed spot-positions are shown in Figure S56. In those Figures, we confirmed that peaks of each resolution distribution are located near zeros, and are formed as nearly Gaussian functions. RMS value represents the positional resolution to 1.5 pixels (100 nm). However, the tail of each distribution appears to be asymmetric. One of the possible explanations of the asymmetry is because of the z-axis. In our analysis, we assumed that spots are characterized as a 2-D Gaussian function, ignoring the z-axis. 3-D Gaussian fitting may be able to resolve the asymmetry of each distribution.
For each cut, we directly compare the simulated spectra to actual spectra in Figures S57 to S60. The number of detected spots per frame per cell for each cuts are shown in Table S14. The left panel shows that the simulated distribution of reconstructed parameters is directly compared with Hiroshima's dataset. The right panel shows the ratio of the simulated spectra to the actual ones in each of bins. The errors are not only observed statistical errors, but also include the simulated sample statistical errors. All simulated spectra are normalized by the number of detected spots.

C.2.2 Experimental configuration A (various ligand inputs)
In this experimental configuration, 10 cell-model samples are prepared for 12 ligand inputs (1.0 pM ∼ 4.0 nM). Cell-compartments for each cell-model are listed in Table S13. We also assume two observational cut-1 and cut-2. Parameter values for the spot-detection algorithm are configured to σ min = 2.0 pixels, σ max = 4.0 pixels, 20 intermediate values in the deviation range, threshold = 15, and overlap = 0.5.
Analysis results are shown in Figures S61; (1) The left and the right top panels show the equilibrium binding curve and the ratio of the reconstructed density to ground-truth. These plots clearly show that shape of the true binding curve is partially recovered in the reconstructed one. While reconstruction efficiency is about 80% and steady at relatively higher ligand concentration inputs (> 0.2 nM), the efficiency is significantly reduced at the lower ligand concentration region (< 0.  Hill coefficient : In a standard approach to biochemistry and biology, the Hill equation (1) can be fitted to the equilibrium binding curves to evaluate signal responses in a receptor system. The minimization function is given by the equation (13). The fitting results are shown in Table S15. While model-truth of the receptor system exhibits the negative cooperativity, coopertivity is positive in the reconstructed curves. Restoration of the cooperative characteristics is thus failed in this measurements.  Table S15: Results of fitting to the Hill equation. The best fit values and uncertainties of each parameters are listed.χ 2 0 is the reduced minimum. If n < 1, then the receptor system increases binding affinity of sites and exhibits negative cooperativity. If n > 1, then cooperativity is positive, decreasing the binding affinity of sites. If n = 1, then there is no cooperativity.

Cut-1
Cut-2 Figure S62: Fitting results. The top panels show the equilibrium binding curves for each cut. The bottom panels show the ratio of the reconstructed binding curves to the fitted curves.
Confidence interval : ∆χ 2 -plots for each area-cuts are shown in Figures S63 show. Black solid lines represent the 1σ (68%), 2σ (95%) and 3σ (99%) confidence intervals. Black and pink points indicate the best fit and ground-truth. The top panels show the ∆χ 2 contour plots: the Hill coefficient (n) vs dissociation constant (K A ). The ground-truth (pink point) is located out of the 3σ confidence contour line, implying the restoration failure of the true parameter values. The middle four plots show the ∆χ 2 plots for each model parameters: the ground-true Hill coefficient and the dissociation constant are clearly out of the contour lines.

Cut-1
Cut-2 Systematic sources : There are three major sources that can generate such systematic shifts. Details are discussed as follows.
(1) Stochasticity: Stochastic process of photon-detection may influence the reconstruction procedures, changing the shape and concavity of the equilibrium binding curve and Scatchard plot. We turn off noise channels in the camera simulation module to generate expected images. Since intensity of the expected images is presented in the level of photoelectron unit, we convert the image intensity unit to the level of ADC counts for proper evaluation. We then applied the analytical procedure to the expected images in order to evaluate the systematic influence on the binding curve and Scatchard plot. First, the expected images are compared with the stochastic images obtained by turning on the noise channels. Figure S64 shows example images of spot-detection on the apical cellsurfaces. The red circles shown in the Figures are more sparsely distributed in the expected images, avoiding capture of the molecular-spots in the defocused region of the expected images. However, more spots are detected near the defocused region of the stochastic images. Figure S65 shows the equilibrium binding curves and Scatchard plot for each area-cuts. The top panels show the comparison of the equilibrium binding curves reconstructed from the expected and stochastic images. The restoration efficiency and defects of the analytical procedure are shown in the middle four panels. We computed the ratio of the reconstructed binding curves of the stochastic images to the expected images. The 2nd panels from the top show the ratio curves, implying that the stochasticity can increase approximately 10% of the restoration efficiency in capturing more molecular-spots. The 2nd panels from the bottom show the fraction of false-spots contaminated in all detected spots, implying that the stochasticity can increase approximately 10% of contamination in capturing the defected spots. Finally, the bottom panels show the comparison of the Scatchard plots reconstructed from the expected and stochastic images. The reconstructed Scatchard plots are aligned in parallel. Consequently, we found no significant change in the shape and concavity of the binding curve and Scatchard plot.
(2) Spot-detection algorithm: Figure S61 shows that the restoration efficiency of the LoG spotdetection method is limited by up to ∼ 80%. The fractional occupancy of the defected spots is steady at ∼ 13% for various ligand concentrations. Such imperfect performance of the spotdetection algorithm may influence to our biological interpretation of cooperativity. Figures S66 and S67 show the time course data in the ligand concentration with inputs more than 0.2 nM. In this concentration range, the reconstructed responses (represented by the red and blue crosses in the Figures) successfully converge to the 80% of the true-equilibrium state (pink dashed lines), reaching to the 80% restoration of the true-equilibrium state. As shown in Figure S61B, the shape and concavity of the reconstructed binding curve and Scatchard plot are unchanged from the ground-truth in the concentration range, implying weak influence on the identification of cooperativity. In addition, more accurate algorithms [19,57,58] may help improving the efficiency and defects, but those algorithms can only help shifting the entire spectra up and down: no significant change in the shape and concavity.
(3) Quasi-static response: The binding system often requires a longer image acquisition period to reach steady state, due to quasi-static response at very low concentration inputs, at ∼ 1 pM scale. This process is a slow transition of the binding system to equilibrium state. If the binding system is quite sensitive to such slow response, then the system cannot converge to the equilibrium state within the acquisition period. The quasi-static responses may accordingly become a critical issue for the biological interpretation of cooperativity. Figures S68 and S69 show the time course data of the ligand concentration with inputs less than 0.2 nM. In this concentration range, the reconstructed responses still remain in a nonequilibrium state within the image acquisition period: 0 to 5, 000 sec. Although the true-nonequilibrium binding state (black lines) is well-restored through the reconstruction procedure, the reconstructed binding state (red and blue crosses) fails to converge to the true-equilibrium states (pink lines), generating a gap between the reconstructed binding state and the true-equilibrium state. Such gaps cannot be even detected during the image acquisition period, thus leading to the misidentification of cooperativity. The ground-true response is represented with a black solid line. The 80% restoration of the true response is shown with a dashed black line.

C.2.3 Experimental configuration A (high concentration range)
In this experimental configuration, we used 10 cell-samples for 6 ligand inputs (200 pM ∼ 4.0 nM). Cell-compartments for each cell-model is listed in Table S13. We also assume two observational cuts: the cut-1 and the cut-2. Parameter values for the spot-detection algorithm are configured to σ min = 2.0 pixels, σ max = 4.0 pixels, 20 intermediate values in the deviation range, threshold = 15, and overlap = 0.5.
Hill coefficient : The Hill equation (1) is fitted to the equilibrium binding curves to evaluate signal responses in the receptor system. The minimization function is given by the equation (13). The fitting results are shown in Table S16 and Figures S70. While model-truth of the receptor system exhibits the negative cooperativity, the Hill coefficients are largely fluctuated around an unity. The cooperative characteristics is thus indeterminable in this measurement.  Table S16: Results of fitting to the Hill equation. The best fit values and uncertainties of each parameters are listed.χ 2 0 is the reduced minimum.
Cut-1 Cut-2 Figure S70: Equilibrium binding curves and Scatchard plot. The top panels show the equilibrium binding curves for each cuts. Bottom ones show the Scatchard plots for each cuts. Blue solid and dashed lines represents the best fit curves and the curves 1σ shifted from the best fits.
Confidence interval : ∆χ 2 -plots for each area-cuts are shown in Figures S71 show. Black solid lines represent the 1σ (68%), 2σ (95%) and 3σ (99%) confidence intervals. Black and pink points indicate the best fit and ground-truth. The top panels show the ∆χ 2 contour plots: the Hill coefficient (n) vs dissociation constant (K A ). The ground-truth (pink point) is located within the 3σ confidence contour line, implying the successful restoration of the true parameter values. The middle four plots show the ∆χ 2 plots for the Hill coefficient and the dissociation constant.

C.2.4 Experimental configuration B (nonequilibrium region)
Single-molecule imaging of apical regions of the dimer model is simulated for the optical specifications and operating conditions of the fluorescence microscopy simulation module shown in Table S6 and S53. Photobleaching is included in the simulated single-molecule images. The results are shown as follows.
(1) Spot-detection : In this experimental configuration, the parameter values for the spot-detection algorithm are configured to σ min = 3.0 pixels, σ max = 6.0 pixels, 30 intermediate values in the deviation range, threshold = 30, and overlap = 0.0. Figure S72 shows example images of spotdetection on the apical cell-regions. Since the focal position is unknown and typically determined by human-eyes, we defined two observational area-cuts (cut-1 and cut-2) with 100 × 300 pixels (132 µm 2 ) represented with dashed lines. While the cut-1 is placed at well-focused region of image center, the cut-2 is shifted out of the image center and located at the defocused region. Red circles represent the spots detected by the LoG method. Each image size is 512 × 512 pixels. Actual minimum and maximum values of the image intensity are 1, 900 and 2, 600 ADC counts. The image intensity is rescaled in the range of 0 to 255.
Analysis results are shown in Table S17 and Figures S73; (1) the Table shows the number and fraction of observed and simulated spots. The simulated spots are true-molecular spots and falsespots that arise from molecules and background noise. In particular, the false-spots are noise-like spots that can mimic the molecular spots. Quite a large fraction of the false-spots are captured with both cuts. (2) The left and the right panels of the Figure show binding curves and fractional occupancy of false-spots for the cut-1 and cut-2. Efficiency of the density reconstruction is approximately 10 % for observational cuts. While the reconstructed area-densities for cuts tend to be linearly increased in time, the observed density is relatively flat in time. In addition, the fractional occupancy of the false-spots is stable in time.   Positional resolution (or localization error) of reconstructed spot-positions are shown in Figure S74. In those Figures, we confirmed that the peaks of each resolution distributions are located near zeros, and are formed as nearly Gaussian distributions. The root mean squared (RMS) value represents the positional resolution to 1.5 pixels (∼ 100 nm). However, the tails of each distributions appear to be asymmetric. One of the possible explanations of the asymmetry relates to the z-axis. In our analysis, we assumed that spots are characterized by a 2-D Gaussian function, ignoring the z-axis. 3-D Gaussian fitting may be able to resolve the asymmetry of each distribution.

Spots/cell Observed Simulated
For each area-cuts, we directly compared the simulated spot-spectra to actual spectra in Figures  S75 and S76. Large discrepancies are clearly observed in all spectra. The left panel shows that the simulated distribution of reconstructed parameters is directly compared with Hiroshima's dataset. The right panel shows the ratio of the simulated spectra to the actual spectra in each of bins. The errors are not only observed statistical errors, but also include the simulated sample statistical errors. All simulated spectra are normalized by the number of detected spots.

Cell-index 1
Hiroshima-2012      x + σ 2 y ) is compared with the actual correlation (right) obtained from Hiroshima's dataset. Comparisons of background pulse-height and SNR distributions are shown in bottom two panels.
(3) Spot-tracking : To avoid linking spots having relatively larger size, we assume maximum threshold in the primary condition of linking two spots. An additional condition is the average size of the two spots must be less than 3.0 pixels, and can be written in the form of σ r i,k + σ r j,k+1 2 < 3.0 pixels (14) where i and j represent spot indexes at k-th and (k + 1)-th image-frames.
In the dimer model, 33 and 12 events per cell were reconstructed with each area-cuts. These events can be categorized into two types; (1) Ground-truth : true-monomer (R, rR, r'R) and -dimer (RR) events are the event-trajectories that are truly reconstructed with true-monomer and -dimer spots. State-transition events are the events that can interact between true-monomer and true-dimer events. While tracking spots, true-monomer (or true-dimer) can often interact with other-type of molecules. For example, the first and second panels of Figure S83 show interactions of molecular-states. Red and blue spots represent true-monomer and true-dimer. Such state-transition events cannot be distinguished from other events. (2) Failure : unfortunately 70-80% of reconstructed event-trajectories are false-events. Event-trajectories are reconstructed under the assumption that the spot-tracking algorithm can successfully track identical molecules (or same molecular ID assigned by Spatiocyte). If the spot-tracking algorithm fails tracking the same molecule and captures other molecules traveling nearby, then the reconstructed event-trajectory is considered to be a false-event. Also, the reconstructed event-trajectory contaminated with non-molecular spots (or false-spots) is consider to be a false-event.
Properties of event-trajectories are shown in Table S18 and Figure S79. The Table shows