Simulation of live-cell imaging system reveals hidden uncertainties in cooperative binding measurements
Masaki Watabe Satya N. V. Arjunan Wei Xiang Chew Kazunari Kaizu Koichi Takahashi

 Cloth Manipulation Using Random-Forest-Based Imitation Learning

    Biao Jia, Zherong Pan, Zhe Hu, Jia Pan, Dinesh Manocha 

<http://cs.unc.edu/ biao/robustm>
Biao Jia is with the Department of Computer Science, University of Maryland at College Park. E-mail:  biao@cs.umd.edu.Zherong Pan is with the Department of Computer Science, University of North Carolina at Chapel Hill. E-mail:  zherong@cs.unc.edu.Zhe Hu is with the Department of Mechanical Biomedical Engineering, City University of Hong Kong.
 E-mail:  zhe.hu@my.cityu.edu.hk.Jia Pan is with the Department of Computer Science, the University of Hong Kong. E-mail:  jpan@cs.hku.hk.Dinesh Manocha is with Departments of Computer Science and Electrical & Computer Engineering, University of Maryland at College Park. E-mail:  dm@cs.umd.edu.
    
===============================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================






We propose a computational method to quantitatively evaluate the systematic uncertainties that arise from undetectable sources in biological measurements using live-cell imaging techniques. We then demonstrate this method in measuring biological cooperativity of molecular binding networks: in particular, ligand molecules binding to cell surface receptor proteins. Our results show how the non-statistical uncertainties lead to invalid identification of the measured cooperativity. Through this computational scheme, the biological interpretation can be more objectively evaluated and understood under a specific experimental configuration of interest.



Keywords : error analysis, systematic uncertainty, bioimage simulation, single-molecule experiment, fluorescence microscopy, cooperative binding, cooperativity



URL: <https://link.aps.org/doi/10.1103/PhysRevE.100.010402> 

DOI: 10.1103/PhysRevE.100.010402 

E-mail: masaki@riken.jp



















§ INTRODUCTION



  

Recent progress in robotic and automated techniques for sampling and analyzing complex biological data can reduce statistical uncertainties, increasing precisions in measuring biological and physical properties in living cells <cit.>. However, despite advances in computational techniques, the measurement process has resisted quantitative evaluation of systematic uncertainties that arise from inaccuracies in experimental data acquisition and analysis <cit.>. An absence and ignorance of evaluation of such systematic uncertainties often lead to excessive interpretations of the measured properties.



  

A key challenge to evaluating the systematic uncertainties is finding meaningful and nonintuitive variance. There exist a large number of systematic sources in the measurement processes, but no well-defined procedure to evaluate systematic variance. Although experienced experimental biophysicists are able to anticipate some systematic sources (e.g., the faulty calibration of measurement equipment) and ensure that most systematic uncertainties are much less than the required precision, other sources cannot even be detected by empirical approaches. The systematic variance that arises from undetectable sources in the measurement process can cause erroneous identification and interpretation of measured properties <cit.>. For example, structural uncertainties that arise from various model assumptions in biological network topologies cannot be directly extracted from experimental data, introducing errors into analysis and potentially giving rise to misleading conclusions <cit.>. Thus, the measurement processes require a non-experimental evaluation method that allows biophysicists to objectively interpret the measured outputs and draw proper conclusions.



  

To better understand the origin of these hidden uncertainties, we consider a process of data-driven (or inductive) modeling in bioimaging. For the sake of simplicity, we assume modeling the steady state behavior of a specific type of protein (e.g., cell-membrane receptors) in a living cell of interest. In this modeling, we also assume that a fluorescence microscopy system measures the concentrations of the proteins tagged with fluorescent emitters (e.g., green fluorescence proteins) within the focal plane of the optics onto a digital camera placed at the conjugated focal plane. Technical details of extracting the protein concentrations from microscopy images are also omitted in this simplification. 








  

Here, as an example, we evaluate the systematic uncertainties associated with a specific parameter (e.g., image acquisition periods) in the simplified measurement process. We first consider that the measured concentration of the observed proteins is a function of time at the i-th image-frame C (t_i). For m cell-samples, the experimenters ensure the stability of the concentration changes in an image acquisition period, T. The observed rate in the protein concentrations converges near 0 as t_i → T,

    Δ C_i/Δ t_i→ 0 ±σ_ C'

where Δ C_i = C (t_i) - C (t_i-1) and Δ t_i = t_i - t_i-1 are the concentration and time difference at the i-th image frame, respectively. σ_ C' represents the statistical deviations in the observed rates. An “apparent" steady state of the protein concentrations is then computed by the time-average integration C±σ_ C.



  

In the data-driven approach, the measured protein concentrations are typically fitted to the network model constructed under the data-interpretation that the measured concentrations fully converge to equilibrium within the acquisition period. The rate of concentration changes in the observed proteins interacting with N variables can be theoretically modeled in the form of an ordinary differential equation. In equilibrium modeling, the rate fully converges to 0 as t →∞,

    d R (t)/dt = f (R(t), x(t),t; θ)→ 0

where R (t) is the observed protein concentration as a function of time and model-parameters θ (e.g., dissociation constants), x (t) = [x_ 1 (t), ... , x_ N (t)] and x_ i (t) is the value of the i-th variable at t and θ. In the equilibrium modeling, the rate of changes in the i-th variable over time also fully converges to 0 as t →∞, d x_i (t)/dt = f_i (R(t), x(t),t; θ)→ 0. The full-equilibrium concentration in the observed proteins is then computed by the time-average integration, R (θ) = lim_τ→∞1/τ∫^τ_0 R (t; θ)dt. Fitting procedures (e.g., maximum-likelihood method) directly compare the measured concentration (C) to the theoretical full-equilibrium concentration (R (θ)) binned into each histogram. The best-fit model parameters (θ) are those which minimize the discrepancies defined in an optimization function ℳ( R (θ), C, σ_ C). Confidence levels and correlations are also estimated in each fitting parameter.



  

Although the experimenters ensure the stability of the measured protein concentration changes in time and show the goodness of fit to the equilibrium model, actual biological cells generally operate out of equilibrium. In Eq. (<ref>), the apparent state stability can be not only interpreted as the complete convergence of the protein responses to full-equilibrium but also the incomplete convergence if the responses are so slow that the protein concentration remains nonequilibrium during the acquisition period. To split this double data-interpretations, the experimenters, however, must contend with the experimental configuration and the state-transition speed in the observed proteins. The data-interpretation that gives rise to the complete convergence requires not only maximizing event samples in a shorter image acquisition period but also fasten the observed rates of converging to the full-equilibrium. The configuration to reduce statistical uncertainties constrains the limitation and sensitivity in measuring the state-transition time. Because of these tradeoffs, it is unclear when or whether the measured protein concentration actually converges to the full-equilibrium within the acquisition period T. 



  

To see the effects of the double data-interpretations in the apparent steady states, we consider a Taylor expansion of Eq. (<ref>) at the model-true equilibrium concentration R = R(θ),

    f (R(t; θ))≈( R (t; θ) - R) ∂ f/∂ R|_ R (t; θ) = R + ⋯

where zeroth-order term vanishes at the model-true equilibrium f (R) = 0. At t = T, we evaluate the systematic variance that arises from the following model assumptions. If full-equilibrium is proper as a model assumption, then f (R (T; θ)) can fully converge to 0. No systematic variance can be generated in the fitting results, | R (T; θ) - R| < σ_C, implying successful restoration of the model-true equilibrium concentrations. However, if nonequilibrium is the proper assumption, then the f (R (T; θ)) can converge to a finite value near 0. Under this assumption, measuring the model-true equilibrium requires a limit excess of the experimental configuration, thereby generating an undetectable gap in the fitting results, | R (T; θ) - R| ≫σ_C. This implies the reconstruction failure of the model-true equilibrium concentration at T, leading to excessive interpretations of the measured protein concentration. Thus, the systematic variance that arises from the double data-interpretations cannot be evaluated without computing the dynamical behavior of the equilibrium models. 



  

A computational modeling approach for a whole experimental system is more relevant to extract the impact of various systematic uncertainties <cit.>. In this work, we introduce a comprehensive method to quantitatively evaluate the systematic variance 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 bioimage simulation module for an oblique illumination fluorescence microscopy system configured to observe receptor proteins binding to ligands on an apical surface of biological cells. We then show how the non-statistical variance leads to misidentification of cooperativity in the binding system.









































































































































































































































§ COMPUTATIONAL METHOD



  

A key insight into evaluating systematic variance is to estimate how well the ground-true model properties can be restored through analytical procedures, influencing the interpretation of biological properties reconstructed (or extracted) from actual biological images <cit.>. Such model-driven evaluation allows us to quantify the restoration efficiency and defects (or failure) in the reconstruction processes. We cannot fully know the true-model; otherwise, there would not be any uncertainty in the biological measurements. However, we can approximately know what to expect either from earlier experimental results or from the biological models derived from experimental knowledge. Such approximations can function as a guide to compute systematic variance in an organized manner from a specific experimental configuration of interest.



  

Systematic variance generally causes the reconstructed properties to be shifted in one direction from the ground-true model properties <cit.>. Of particular importance is the quantifying of systematic variance that arises from undetectable sources in the measurement processes. If the ground-true property is well-restored through the reconstruction process, then it is unchanged and has weak influence on the biological interpretation. However, if the reconstruction process poorly restores the true property, it may change significantly, affecting the biological interpretation. For example, geometric uncertainties imposed by the unobservable sources such as complex and irregular shapes of membrane-enclosed cellular compartments (e.g., endoplasmic reticulum) in FRAP experiments lead to invalid measurement of molecular mobility and compartmental connectivity <cit.>.



  

A realistic ground-true model simulation that approximately represents a whole experimental system can be constructed for the model-driven evaluation. Biological measurements using live-cell imaging techniques are generally governed by various natural laws and principles of biochemistry and physics. Models of each imaging process can be simulated within the limited range and dimensions of model-parameter space. Various model-simulation studies exist in physics simulations for molecular fluorescence and optical apparatus <cit.> and systems biology simulations aimed at explaining and predicting biological phenomenon <cit.>. We integrate these model-simulations into a unified model corresponding approximately to a whole experimental system, thus helping to explore the extended dimensions of model-parameter space that can affect potential imaging and analytical outcomes <cit.>. In particular, we have developed the bioimage simulation platform for handling a large range of biochemical and physical parameters that govern image-based measurement systems, generating computational photomicrographs that arise from the various systematic sources: spatiotemporal model of biological cells, photophysics and imaging apparatus <cit.>. Through such simulation platforms, the biological interpretation of measured properties 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 single-molecule images in an equilibrium region, and then identified apparent (or observed) cooperativity to be "negative" for the dimer formation <cit.>. 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 <cit.>. 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 <cit.>. In particular, negative cooperativity has been incorporated in the model by introducing an intermediate state transition. Figure <ref>A illustrates the biochemical reaction network of the dimer formation. Values of model parameters are shown in Table <ref>. 

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 <ref> illustrates the biochemical reaction network of the higher-order oligomer formation. Values of additional model parameters are shown in Table <ref>.




  

Figures <ref>B 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 <ref>. 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 <ref>. 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 <cit.> 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 <ref>A and B show the comparison of ground-true spot area-density (2.372  spots/μ m^2) to the reconstructed one. Figure <ref>A 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 <ref>B 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 <ref>C 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 <ref>A shows the difference between the ground-true spot-position (x_0^ true, y_0^ true) and the reconstructed one (x_0^ reco, y_0^ reco). In the 2-dimensional distribution of the localization errors: (x_0^ reco - x_0^ true, y_0^ reco - y_0^ 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 <ref>B 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: light-scattering, autofluorescence, and optical aberration. More detailed implementation of the simulation module is required to resolve the discrepancies. Furthermore, such direct comparisons for the spot-properties 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 <ref>A 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). Figure<ref>C 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

    B(L)     =    B_0 L^n/K^n_A + L^n

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 <ref>. 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 <ref>D shows the comparison of the true Scatchard plot to the reconstructed one. The plot shows that the concavity of the true plot is not 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 <ref>A 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 <ref>C 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 true-equilibrium state (pink dashed lines), reaching the 80% restoration of the true-equilibrium state. As we showed in Figure <ref>B, 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 <ref>E 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 <cit.>. Figure <ref>A 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 <ref> and in Figures <ref>. Figure <ref>B 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



  

Major scientific activities in biological sciences are dedicated not only to extracting laws and patterns from experimental data but also experimentally validating the biological models derived from experimental knowledge <cit.>. Various model candidates can be either confirmed or refuted by repeating these activities. The key to successful reporting of experimental results is to provide an objective evaluation and representation of the uncertainties that arise from imprecision and inaccuracies in the experimental processes. The study and estimation of the experimental uncertainties have been generally known as error analysis, its main function being to allow biophysicists to numerically indicate the validity and confidence of their experimental results <cit.>.



  

In the error analysis, statistical analysis of the experimental data is only half the story. The other half is the computation of the systematic uncertainties that affect the sensitivity and limitation of a given experimental configuration to the model parameters and, even more significantly, to new advancements in biophysics and biology. An identification of the new findings in the experimental approaches must always contend with those error estimations. For this reason, we proposed the computational method to evaluate the impact of various systematic uncertainties in the biological measurements using live-cell imaging techniques. We then presented the first examples of not only estimating the systematic uncertainties in model assumptions and parameters (e.g., image acquisition periods and spot-detection efficiency) that can affect the cooperative binding measurements but also reducing them to levels allowing for proper conclusions. In the near future, we believe our computational scheme can help bridges the gap between theory and experiment in biological sciences.














§ ACKNOWLEDGEMENTS

We would like to thank Yasushi Okada, Tomonobu M. Watanabe, Jun Kozuka, Michio Hiroshima, Kozo Nishida, Kazunari Iwamoto, Yuki Shindo, Hanae Shimo, Yin Fai Chin, Suguru Kato, Toru Niina, Shin Moriga, Taku Tsuzuki, Koji Ochiai, Keiko Itano, Sibel Ozer, Kotone Itaya and Kaoru Ikegami for their guidance and support throughout this research work. We would also like to thank Yasushi Sako and Kylius Wilkins for their critical reading of the manuscript. This research work is supported by the JSPS (Japanese Society for the Promotion of Science) KAKENHI Grants-in-Aid program for Challenging Exploratory Research (15K12146).



§ AUTHOR CONTRIBUTIONS

MW conceived and designed the problem. MW programmed the fluorescence microscopy simulation module. SNVA, WXC and KK develop E-Cell simulator. MW constructed cell-models. MW performed analyses. MW wrote manuscript. KK and KT provided support and guidance.
























































































































10

king2009
King RD et al. (2009) The automation of science.
Science (New York, N.Y.) 324(5923):85–89.

marx2013
Marx V (2013) THE BIG CHALLENGES OF BIG DATA.
Nature 498(7453):255–260.

jordan2015
Jordan MI, Mitchell TM (2015) Machine learning: Trends, perspectives, and
  prospects.
Science 349(6245).

yachie2017
Yachie N, Robotic Biology Consortium, Natsume T (2017) Robotic crowd
  biology with Maholo LabDroids.
Nature 35(4):3–5.

amrhein2019
Amrhein V, Mcshane B (2019) Retire statistical significance.
Nature 567(7748):305–307.

geris2016
Geris L, Gemez-Cabrero D (2016) Uncertainty in Biology.
(Springer).

taylor1997
Taylor JR (1997) An Introduction to Error Analysis.
(University Science Book), 2nd edition.

bevington2003
Bevington PR, Robinson DK (2003) Data Reduction and Error Analysis.
(McGrawHill), 3rd edition.

babtie2014
Babtie AC, Kirk P, Stumpf MPH (2014) Topological sensitivity analysis for
  systems biology.
Proc. Natl. Acad. Sci. U.S.A. 111(52):18507–18512.

kirk2015
Kirk PDW, Babtie AC, Stumpf MPH (2015) Systems biology (un)certainties.
Science 350(6259):386–388.

incerti2009
Incerti S et al. (2010) The Geant4-DNA project.
Int. J. Model. Simul. Sci. Comput. 1(2):157–178.

incerti2018
Incerti S,  et al. (2018) Geant4-DNA example applications for track structure
  simulations in liquid water : A report from the Geant4-DNA Project.
Med. Phys. 45(8):e722–e739.

agostinelli2003
Agostinelli S, Allison J, Amako K (2003) Geant4–a simulation toolkit.
Nucl. Instrum. Methods Phys. Res. A 506(506):250–303.

allison2006
Allison J,  et al. (2006) Geant4 developments and applications.
IEEE Trans. Nucl. Sci. 53(1):270–278.

allison2016
Allison J,  et al. (2016) Recent developments in GEANT4.
Nucl. Instrum. Methods Phys. Res. A 835:186–225.

sbalzarini2013
Sbalzarini IF (2013) Modeling and simulation of biological systems from image
  data.
BioEssays : news and reviews in molecular, cellular and
  developmental biology 35(5):482–90.

mai2013
Mai J, Trump S, Lehmann I, Attinger S (2013) Parameter Importance in FRAP
  Acquisition and Analysis: A Simulation Approach.
Biophys. J. 104(9):2089–2097.

lakowicz2006
Lakowicz J (2006) Principles of Fluorescence Spectroscopy.
(Springer), p. 954.

pawley2008
Pawley J (2008) Handbook of Biological Confocal Microscopy.
(Springer), p. 988.

mansuripur2009
Mansuripur M (2009) Classical Optics and its Application.
(Cambridge University Press), p. 714.

valeur2012
Valeur B, Berberan-Santos M (2012) Molecular Fluorescence.
(Wiley), p. 592.

kaneko2006
Kaneko K (2006) Life: An introduction to complex systems biology.
(Springer).

alon2007a
Alon U (2007) An introduction to systems biology — Design principles of
  biological circuits.
(Champman & Hall/CRC).

tomita1999
Tomita M et al. (1999) E-CELL: Software environment for whole-cell
  simulation.
Bioinformatics 15:72–84.

boulanger2009
Boulanger J, Kervrann C, Bouthemy P (2009) A simulation and estimation
  framework for intracellular dynamics and trafficking in video-microscopy and
  fluorescence imagery.
Med. Image Anal. 13(1):132–142.

rezatofighi2013
Rezatofighi SH et al. (2013) A framework for generating realistic synthetic
  sequences of total internal reflection fluorescence microscopy images.
Proceedings - International Symposium on Biomedical Imaging pp.
  157–160.

angiolini2015
Angiolini J, Plachta N, Mocskos E, Levi V (2015) Exploring the Dynamics of
  Cell Processes through Simulations of Fluorescence Microscopy Experiments.
Biophys. J. 108(11):2613–2618.

watabe2015
Watabe M et al. (2015) A Computational Framework for Bioimaging Simulation.
PLOS One 10(7):e0130089.

venkataramani2016
Venkataramani V, Herrmannsdörfer F, Heilemann M, Kuner T (2016) SuReSim:
  simulating localization microscopy experiments from ground truth models.
Nat. Methods 13(4).

linden2016
Lindén M, Ćurić V, Boucharin A, Fange D, Elf J (2016)
  Simulated single molecule microscopy with SMeagol.
Bioinformatics p. btw109.

weigert2018
Weigert M, Subramanian K, Bundschuh ST, Myers W, Kreysing M (2018) Biobeam —
  Multiplexed wave-optical simulations of light-sheet microscopy.
PLOS Comput. Biol. 14(4):e1006079.

hiroshima2012
Hiroshima M, Saeki Y, Okada-Hatakeyama M, Sako Y (2012) Dynamically varying
  interactions between heregulin and ErbB proteins detected by single-molecule
  analysis in living cells.
Proc. Natl. Acad. Sci. U.S.A. 109(35):13984–9.

arjunan2010
Arjunan SNV, Tomita M (2010) A new multicompartmental reaction-diffusion
  modeling method links transient membrane attachment of E. coli MinE to E-ring
  formation.
Syst. Synth. Biol. 4:35–53.

takahashi2005
Takahashi K, Arjunan S, Tomita M (2005) Space in systems biology of signaling
  pathways – towards intracellular molecular crowding in silico.
FEBS letters 579:1783–1788.

vanderwalt2014
van der Walt S et al. (2014) scikit-image: image processing in Python.
PeerJ 2:e453.

kitano2002_sci
Kitano H (2002) Systems biology: a brief overview.
Science (New York, N.Y.) 295(5560):1662–4.

kitano2002_nat
Kitano H (2002) Computational systems biology.
Nature 420(6912):206–10.

godfrey2003
Godfrey-Smith P (2003) Theory and Reality — an introduction to the
  philosophy of science.
(The University of Chicago Press).

zhao2011
Zhao Q, Young IT, de Jong JGS (2011) Photon budget analysis for fluorescence
  lifetime imaging microscopy.
Journal of biomedical optics 16(8):086007.

thompson2010
Thompson Ma, Biteen JS, Lord SJ, Conley NR, Moerner WE (2010) Molecules
  and Methods for Super-Resolution Imaging.
(Elsevier Inc.) Vol. 475, 1 edition, pp. 27–59.

didier2005
Didier P, Guidoni L, Bardou F (2005) Infinite Average Lifetime of an Unstable
  Bright State in the Green Fluorescent Protein.
Physical Review Letters 95(9):090602.

zondervan2004
Zondervan R, Kulzer F, Kol'chenk Ma, Orrit M (2004) Photobleaching of
  Rhodamine 6G in Poly(vinyl alcohol) at the Ensemble and Single-Molecule
  Levels.
The Journal of Physical Chemistry A 108(10):1657–1665.

hoogenboom2007
Hoogenboom JP, Hernando J, Van Dijk EMHP, Van Hulst NF,
  García-Parajó MF (2007) Power-law blinking in the fluorescence
  of single organic molecules.
ChemPhysChem 8(6):823–833.

margolin2005
Margolin G, Barkai E (2005) Nonergodicity of blinking nanocrystals and other
  Lévy-walk processes.
Physical Review Letters 94(8):1–4.

brokmann2003
Brokmann X et al. (2003) Statistical Aging and Nonergodicity in the
  Fluorescence of Single Nanocrystals.
Physical Review Letters 90(12):120601.

shimizu2001
Shimizu K et al. (2001) Blinking statistics in single semiconductor
  nanocrystal quantum dots.
Physical Review B 63(20):205316.

miyanaga2009
Miyanaga Y, Matsuoka S, Ueda M (2009) Single-Molecule Imaging Techniques to
  Visualize Chemotactic Signaling Events on the Membrane of Living
  Dictyostelium Cells.
Methods in Molecular Biology 571:417–435.

axelrod2008
Axelrod D (2008) Chapter 7: Total internal reflection fluorescence
  microscopy.
(Elsevier Inc.) Vol. 89, 1 edition, pp. 169–221.

wazawa2005
Wazawa T, Ueda M, Rietdorf J (2005) Microscopy Techniques - Advances in
  Biochemical Engineering/Biotechnology.
(Springer Berlin / Heidelberg) Vol. 95, pp. 1297–1300.

axelrod2003
Axelrod D (2003) Total internal reflection fluorescence microscopy in cell
  biology.
Methods in enzymology 361(2):1–33.

teramura2006
Teramura Y et al. (2006) Single-molecule analysis of epidermal growth factor
  binding on the surface of living cells.
EMBO J. 25(18):4215–4222.

uyemura2005
Uyemura T, Takagi H, Yanagida T, Sako Y (2005) Single-molecule analysis of
  epidermal growth factor signaling that leads to ultrasensitive calcium
  response.
Biophys. J. 88(5):3720–30.

sako2002
Sako Y, Uyemura T (2002) Total Internal Reflection Fluorescence Microscopy for
  Single-molecule Imaging in Living Cells.
Cell Structure and Function 27(5):357–365.

sako2000
Sako Y, Minoghchi S, Yanagida T (2000) Single-molecule imaging of EGFR
  signalling on the surface of living cells.
Nature Cell Biology 2(3):168–172.

tokunaga1997
Tokunaga M, Kitamura K, Saito K, Iwane AH, Yanagida T (1997) Single molecule
  imaging of fluorophores and enzymatic reactions achieved by objective-type
  total internal reflection fluorescence microscopy.
Biochemical and biophysical research communications
  235(1):47–53.

jena2008
Jena BP (2008) Methods in Nano Cell Biology.
(Academic Press), 1 edition.

basset2015
Basset A, Boulanger J, Salamero J, Bouthemy P, Kervrann C (2015) Adaptive Spot
  Detection With Optimal Scale Selection in Fluorescence Microscopy Images.
IEEE Trans. Image Process. 24(11):4512–4527.

ruusuvuori2010
Ruusuvuori P et al. (2010) Evaluation of methods for detection of fluorescence
  labeled subcellular objects in microscope images.
BMC Bioinform. 11:248.

smal2010
Smal I, Loog M, Niessen W, Meijering E (2010) Quantitative comparison of spot
  detection methods in fluorescence microscopy.
IEEE Trans. Med. Imaging 29(2):282–301.







 
 

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.






toc.2pt





§ IMPLEMENTATION UPDATES


r9.3cm


  
    < g r a p h i c s >

  [Optical configurations]Optical configurations of the TIRFM simulation module <cit.>.
  






  

We have implemented the simulation modules of total internal reflection fluorescence microscopy (TIRFM) and laser-scanning confocal microscopy (LSCM) <cit.>. 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. Epi-illumination and oblique illumination techniques are implemented. 





 §.§ 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 <ref> <cit.>. Implementation assumptions are summarized in Table <ref>. 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. 










  §.§.§ Molecular fluorescence



  

Molecular fluorescence is the result of physical and chemical processes in which fluorophores emit photons from electronically excited states <cit.>. 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 <cit.> : 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 <ref> 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 

    n_1    =     n_0× 10^-A
           where   n_0 = σ δ T/4 π|  A_T|^2

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)/N_A ϵ 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

    A =  log( I/I_0) = ϵ c l

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

    n_emit    =    Φ_QY n_abs
        =    Φ_QY( n_0 - n_1) 
        =    Φ_QY n_0( 1 - 10^-A)

In previous implementations of bioimaging simulation <cit.>, 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

    n_emit   ≈   Φ_QY n_0 ln(10) A 
       ∼   Φ_QY n_0× 10^-6

The right panel of Figure <ref> represents the number of emitted photons as a function of absorption coefficient for various flux densities. Red and black lines represent Beer-Lamberts law and its approximation for various incident beam flux densities. The intersection of the black and red lines is found in the typical range of the absorption coefficients; 10^4 to 10^6  [1/(M cm)]. 


  (2) Aging Mechanism <cit.> : 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

    P(τ) = α/τ_0( τ_0/τ)^1 + α

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.




  (3) Photoblinking <cit.> : Time series of each fluorophore often exhibit fluorescence intermittency or blinking, where at random times the fluorophores display the transition between a bright state (ON state in which it stays a time τ_on) and a reversible dark state (OFF state in which it stays a time τ_off). The process is characterized based on a time sequence of ON and OFF {τ^(1)_on, τ^(1)_off, τ^(2)_ on, τ^(2)_ off, τ^(3)_on, τ^(3)_off, ..... τ^(n)_on, τ^(n)_off}.  The time τ^(i)_on/off are drawn at random from the probability distribution given by the equation (<ref>).










  §.§.§ Examples of single-molecule images



  

We constructed realistic particle model of TMR-tagged proteins on glass surface as shown in Figure <ref>. 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 <ref>. Incident beam angle is 65.7^∘ (Evanescent fields). Photophysical parameterizations of fluorophores are also shown in the Table. Results are shown in Figures <ref> and <ref>. These Figures show snapshot images captured with TIRFM imaging configuration.





















 §.§ 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 <ref>. Implementation assumptions are summarized in Table <ref>. 





  §.§.§ 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 <ref> illustrates schematic views of epi-illumination, oblique illumination and Evanescent field. The left panel of Figure <ref> 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 <ref> shows the penetration depth of the evanescent field as a function of the incident beam angle. More details are described in reference <cit.>. 





  
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 <ref>. Figure <ref> 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 <ref>-<ref>. Each figures show snapshot images captured with fluorescence miroscopy imaging configuration.















  §.§.§ 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 <cit.>. This illumination technique allows us to investigate the dynamics of individual molecules at an apical cell surface. The left panel of Figure <ref> 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) <cit.>, 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 <ref>. 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 <ref> shows relative intensity of incident electric fields as a function of incident beam angle. Results are shown in Figures <ref> and <ref>. 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 <cit.>. Figure <ref> 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 <ref>. Results are shown in Figures <ref> and <ref>. Each Figure shows snapshot images changing focal heights from top regions of the cell-model. 





























§ MODEL CONSTRUCTION



 §.§ 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 <cit.>. The left panel of Figure <ref> 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 <ref>. We simulated the dimerization model without the systematic effects that arise from microscopy specification and its operating condition. Simulation results are shown in Figure <ref>. 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 <cit.>. 





 §.§ 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 <ref>B illustrates an additional reaction network of the higher-order oligomer formation. Values of each parameter are shown in Table <ref>. 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) <cit.>. 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). 









  
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 <ref>. Comparisons of actual single-molecule images to the simulated single-molecule images for 1.0  nM ligand input are shown in Figure <ref>-<ref>. The left column shows the actual images obtained by the single-molecule experiment (Movie S1 in Ref. <cit.>). 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. 



















 §.§ PSF-verification



  

Verification is the process of confirming the simulation modules are correctly implemented with respect to conceptual description and analytical solutions <cit.>. 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 <ref>. 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 <ref> show the comparison of the simulated-form (red lines) to 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 simulated-form of the PSF. The top panels of Figures <ref> 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 <ref> and <ref> 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.



































































§ 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 <cit.>[Also see <https://en.wikipedia.org/wiki/Blob_detection> and <http://scikit-image.org/docs/dev/auto_examples/features_detection/plot_blob.html>.]. The method detects spot-like features by searching for scale space extrema of a scale-normalized LoG :

    ∇^2_norm L(x, y; σ)     =    σ^2 ( L_xx + L_yy)

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. <cit.>.

  (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

    G(x,y)     =     N_0 exp( - (x - x_0)^2/2 σ^2_x - (y - y_0)^2/2 σ^2_y) + bg

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

    | r⃗_j,k+1 - r⃗_i,k|     <    σ_r_i,k + σ_r_j,k+1/2

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. 







 §.§ Simple cell-models



  §.§.§ 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 <ref>. Photobleaching is included in the simulated single-molecule 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 <ref>.





  
(1) Spot-detection :
In this test 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. Figure <ref> 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 <ref> shows the reconstruction efficiency of the area-density in apical and basal cell-regions. Ground-true spot area-density 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.





  
(2) Spot-property :
Using the simulated images including photobleaching effects, we reconstruct spot-properties. The reconstructed spot-properties are presented as six parameters of the Gaussian function; spot pulse-height (or normalization factor), central position, spot pulse-width and background pulse-height. Distributions and correlations of each parameters are shown in Figures <ref>, <ref> and <ref>.


  (a) Spot pulse-height (N_0) : Distributions of the spot pulse-height are shown in the top panels of Figure <ref>. This is also known as cluster size in Hiroshima's work <cit.>. 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 monomer-like spots. The histograms show that the photobleached-dimer spots (blue) are identified as monomer-like spots (red).

  (b) Localization errors (r⃗^ reco - r⃗^ true) : Positional resolutions (or localization error) of reconstructed spot-positions are shown in the bottom three panels of Figure <ref>. 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.

  (c) Spot-size (σ_x, σ_y and σ_x vs σ_y) : Distributions of spot-size (or Gaussian-widths) are shown in top three panels of Figure <ref>. Those histograms show that the spots are nearly symmetric in x-y axes and their-size is about 1.35  pixels (89.63  nm).

  (d) Background pulse-height (bg) : Distribution of background pulse-height is shown in the bottom panels of Figures <ref>. The background pulse-heights are distributed around the offset value, and nearly form a Gaussian distribution. 

  (e) SNR : The top panels of Figures <ref> show SNR distributions. We simply assume that the SNR is the ratio of the spot pulse-height to the spot and background noise. It can be written in the form of SNR = s/√(s + b) where s and b are spot and background pulse-height in photoelectron unit. The overall SNR distribution of the basal cell-images is shifted to relatively larger SNR values than that of the apical cell-images. 

  (f) Correlations : Correlation distributions of the spot pulse-height to spot-size and localization error are shown in the middle and the bottom panels of Figures <ref>. Two peaks are observed in those correlations. First and second peaks correspond to monomer-like and dimer-like spots. 
























  
(3) Spot-tracking :
In the simple cell-model, 118 and 382 events are reconstructed in the apical and basal cell-surfaces. These events can be categorized into three types; true-monomer, true-dimer, and false event-reconstruction. Table <ref> shows the number and fraction of each event types. Figures <ref> show space-time series of the observed and ground-truth reconstructed event-trajectories in the apical and basal cell-surfaces. Black spots represent observed states of detected spots. Red and blue spots represent the ground-truth of monomer and dimer. 

Figures <ref> represent all reconstructed trajectories for each event-types. Figures <ref> to <ref> show example trajectories of each event-types and their reconstructed intensity (or spot pulse-height) changes in time; (a) True-monomer events are successfully reconstructed with lower-intensity. (b) True-dimer events are successfully reconstructed with higher-intensity. (c) Event-reconstruction failed. While tracking spots, true-monomer (or true-dimer) can be often exchanged with a neighbour molecule. For example, the third and fourth panels of Figure <ref> and <ref> show exchanges of molecular-states and molecular identification-numbers assigned by Spatiocyte cell simulations. Such event-exchanges cannot be distinguished from other events. 





























  
(4) Event-property :
Reconstructed track properties are presented as physical parameters: event-length, event-vertex, diffusion-coefficient, and event pulse-height. Distributions of each physical parameters are shown in Figures <ref>, <ref> and <ref>. 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 <cit.>. 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

    < Δ r^2(δ t) >     =     4 D δ t + a

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 <ref>. 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 <ref>. 

  (c) Event pulse-height : Distributions of event pulse-heights are shown in the top panels of Figure <ref>. The event pulse-height is defined as the total spot pulse-height in a given trajectory, and written in the form of

    N_evt    =    ∑^spots_j = 0 N_0 (r⃗_j)

where N_0 is the spot pulse-height of j-th spot in a given trajectory. In addition, middle two panels of the Figures show distributions of event pulse-height per path-length and time-frame. Separation of true-monomer and true-dimer event ditributions can be clearly seen in those Figures.

  (d) Event-vertex : Initial and final event-vertex distributions are shown in Figure <ref>. Displacement distribution of the initial vertex to the final vertex is also shown in the bottom panels of Figure <ref>. Initial and final x- and y-vertex distributions are expected to be flat in the apical and basal cell-regions. However, the vertex distributions are not flat in those Figures. More cell-model samples are required to verify the uniformity of each vertex distributions.


















  §.§.§ 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 <ref>. Experimental configuration is shown in Table <ref>. 



Experimental configuration A (1.0  nM ligand input)


  

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 <ref>. 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 <ref> show example images of spot-detection 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 <ref> and Figures <ref>; (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 Gaussian-width) and background pulse-height. Distributions and correlations for each parameter is shown in Figures <ref> to <ref>. 






















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 <ref>. 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 <ref>. Red and blue crosses represent the reconstructed values obtained by the cut-1 and cut-2. Ground-truth and its 50% recovery are presented with black solid and dashed lines. (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

    χ^2     =    ∑^N bins_i = 0(E_i - O_i)^2/σ_i^2

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 <ref>. 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. 











  
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 <ref> 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 <ref> 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.





r10.3cm


  
    < g r a p h i c s >

    
  Experimental configurations <cit.>.
  






 §.§ 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 <ref> and <ref>. 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 <cit.>. Details of each experimental configurations are shown in Figure <ref> and Table <ref>. 








  §.§.§ 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 <ref>. In particular, photobleaching is not included in this experimental configuration. 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 = 2.0  pixels, σ_max = 4.0  pixels, 20 intermediate values in the deviation range, threshold = 15, and overlap = 0.5. Figures <ref> show example images of spot-detection on the apical cell-regions. We assume two observational area-cuts (cut-1 and cut-2) to 80 × 80  pixels (28 μ m^2) and 200 × 200  pixels (176 μ m^2) represented with inner and outer dashed boxes. 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,500 ADC counts. The image intensity is rescaled in the range of 0 to 255. The microscopy is configured to focus on the image center. While molecular-spots are clear and properly picked up at the focused area, the spots are blured and not well detected near the curved-region of cellular membrane. 



  

Analysis results of various observational area-cuts are shown in Table <ref> and Figures <ref>; (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. 


















  
(2) Spot-property :
Reconstructed spot-properties are presented as the Gaussian function of six parameters; spot pulse-height (or normalization factor), central position, spot-size (or Gaussian-width) and background pulse-height. Distributions and correlations of each parameter are shown in Figures <ref> to <ref>.



  

Positional resolutions (or localization error) of reconstructed spot-positions are shown in Figure <ref>. 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 <ref> to <ref>. The number of detected spots per frame per cell for each cuts are shown in Table <ref>. 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.



























  §.§.§ 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 <ref>. 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 <ref>; (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.2  nM). (2) The bottom left panel show fraction of false-spots contaminating all detected spots. Approximately 13% of the detected spots are false-spots with the cut-1. Roughly, 70% 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 concave-up curve. However, shape of the reconstructed plots is a concave-down curve or a straight line. Thus, measurement effect can significantly change apparent biological properties in single-molecule experiments.







  
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 <ref>. 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.















  
Confidence interval :
Δχ^2-plots for each area-cuts are shown in Figures <ref> 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. 







  
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 <ref> shows example images of spot-detection on the apical cell-surfaces. 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 <ref> 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 <ref> shows that the restoration efficiency of the LoG spot-detection 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 spot-detection algorithm may influence to our biological interpretation of cooperativity.

Figures <ref> and <ref> 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 <ref>B, 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 <cit.> 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 <ref> and <ref> 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.


























  §.§.§ 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 <ref>. 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 <ref> and Figures <ref>. 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.

















  
Confidence interval :
Δχ^2-plots for each area-cuts are shown in Figures <ref> 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.







  §.§.§ 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 <ref> and <ref>. 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 <ref> shows example images of spot-detection 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 <ref> and Figures <ref>; (1) 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. 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. 





  
(2) Spot-property :
Reconstructed spot-properties are presented as the Gaussian function of six parameters; spot pulse-height (or normalization factor), central position, spot-size (or Gaussian-width) and background pulse-height. Distributions and correlations of each parameter are shown in Figures <ref> to <ref>. 



  

Positional resolution (or localization error) of reconstructed spot-positions are shown in Figure <ref>. 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.



  

For each area-cuts, we directly compared the simulated spot-spectra to actual spectra in Figures <ref> and <ref>. 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. 



























  
(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

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 <ref> 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 <ref> and Figure <ref>. The Table shows the number and fraction of each event types. The Figures show space-time series of the observed and simulated event-trajectories. Black and red spots represent the spots selected by the cuts and example event-trajectories. Figures <ref> to <ref> show examples of the observed and simulated event-trajectories and their reconstructed intensity (or spot pulse-height) changes in time. 





























  
(4) Event-property :
For each area-cut, we directly compared the simulated event-spectra to actual spectra in Figures <ref> to <ref>. Red, blue, yellow and green filled histograms represent event-distribution of true-monomer, true-dimer, state-transition and false-reconstruction events. Black solid lines represent the sum of true-distributions. Black cross represents actual event-distributions. The left panel shows  the simulated distribution of reconstructed parameters directly compared with Hiroshima's dataset. The right panel shows the ratio of the simulated spectra to the actual spectra in each bin. 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. 















