## New Research In

### Physical Sciences

### Social Sciences

#### Featured Portals

#### Articles by Topic

### Biological Sciences

#### Featured Portals

#### Articles by Topic

- Agricultural Sciences
- Anthropology
- Applied Biological Sciences
- Biochemistry
- Biophysics and Computational Biology
- Cell Biology
- Developmental Biology
- Ecology
- Environmental Sciences
- Evolution
- Genetics
- Immunology and Inflammation
- Medical Sciences
- Microbiology
- Neuroscience
- Pharmacology
- Physiology
- Plant Biology
- Population Biology
- Psychological and Cognitive Sciences
- Sustainability Science
- Systems Biology

# Stochastic approach to the molecular counting problem in superresolution microscopy

Edited by R. Stephen Berry, University of Chicago, Chicago, IL, and approved November 19, 2014 (received for review May 1, 2014)

## Significance

Large and complex macromolecular assemblies—like RNA and DNA polymerases, the kinetochore, the ribosome, ATP synthase, and many others—are critical components of the cell. To fully characterize these molecular machines, it is not sufficient to rely solely on traditional structural methods, like cryo-EM and X-ray crystallography that provide great structural detail in vitro; we need experimental and theoretical tools that can describe the organization of these machines in their native cellular environment so as to better understand their function. Here we provide a strategy for extracting precisely this information directly from superresolution imaging data, a state-of-the-art technique for probing biological structures in living cells below the diffraction limit.

## Abstract

Superresolution imaging methods—now widely used to characterize biological structures below the diffraction limit—are poised to reveal in quantitative detail the stoichiometry of protein complexes in living cells. In practice, the photophysical properties of the fluorophores used as tags in superresolution methods have posed a severe theoretical challenge toward achieving this goal. Here we develop a stochastic approach to enumerate fluorophores in a diffraction-limited area measured by superresolution microscopy. The method is a generalization of aggregated Markov methods developed in the ion channel literature for studying gating dynamics. We show that the method accurately and precisely enumerates fluorophores in simulated data while simultaneously determining the kinetic rates that govern the stochastic photophysics of the fluorophores to improve the prediction’s accuracy. This stochastic method overcomes several critical limitations of temporal thresholding methods.

Protein–protein and protein–nucleic acid interactions are responsible for most of information processing and control in the cell. Moreover, essential cellular tasks such as replication, transcription, translation, recombination, control of gene expression, and transport of proteins across membranes, to cite a few, all depend on the interaction among preorganized molecular assemblies (1⇓⇓–4). Developing a mechanistic understanding of cell biology requires a quantitative determination of macromolecular organization and interaction in living cells. In particular, characterizing protein–protein and protein–nucleic acid complexes and their stoichiometric ratios is a first step in determining how altered stoichiometries may lead to disease states of the cell (5⇓⇓–8).

Conventional optical microscopy is diffraction limited and typically cannot resolve images below the 250-nm range, but superresolution (SR) methods can achieve tens of nanometer resolution (9⇓⇓–12). One such SR method is photoactivated localization microscopy (PALM), where temporal separation is achieved by illuminating a sample with inactive (i.e., nonfluorescing) fluorophores under low light intensity. The light stochastically triggers fluorophore activation. Once active, a fluorophore is excited by light of a different wavelength and releases a burst of photons (an emission burst). This emission is used to locate the position of the emitter with uncertainty that depends only on the number of photons detected. A short time later, the fluorophore irreversibly photobleaches. The low probability of photoactivation ensures that two fluorophores separated less than the diffraction limit will most probably not emit simultaneously. In this manner, the PALM fluorophores that could not otherwise be spatially separated are instead separated in time. The fluorophores in PALM are genetically encoded with photoactivatable fluorescent proteins (PA-FPs) (13)—often mEos2 (14) or Dendra2 (15, 16)—that are fused to proteins of interest.

Generating PALM images is now routine, and both SR tagging techniques and instrumentation have progressed much faster than our ability to analyze PALM images. In particular, extracting protein complex stoichiometry information from PALM data is not straightforward. Rather, stoichiometry is often determined using in vitro (e.g., coimmunoprecipitation) (17) or coarser (diffraction-limited) methods (18⇓⇓⇓⇓–23). However, PALM has the potential to provide molecular counting with single molecule sensitivity and in vivo. However, several obstacles remain: (*i*) PA-FP “blinking” leads to severe overcounting biases. Blinking refers to a process by which a PA-FP produces a series of intermittent emission bursts, instead of a single continuous burst (24, 25), by transiently transitioning to a “dark” state. Accordingly, counting the number of emission bursts over a region of interest (ROI) will lead to overestimating the number of labeled molecules; and (*ii*) Unknown blinking statistics. Often the blinking properties of common PA-FPs are known in vitro, but not in vivo in the specific cellular compartment occupied by the molecules. Characterizing the blinking properties of PA-FPs in their cellular context is required because these properties are highly sensitive to their local cellular environment (pH, ionic strength) (26⇓⇓⇓⇓⇓⇓–33). An analytical method capable of extracting the blinking statistics of PA-FPs inside the cell to correct for blinking is not available.

Here, we describe a stochastic approach for dealing with these obstacles. Our approach is an adaptation of the continuous time aggregated Markov model (AMM) techniques developed in the ion channel literature to estimate kinetic rates using maximum-likelihood methods for channel opening and closing events in patch-clamp experiments (34⇓⇓⇓–38). Our objective is, thus, to arrive at molecular quantification in PALM. Previous studies have addressed the PALM counting problem by setting a temporal threshold

In what follows we develop the theoretical methodology and benchmark the method on sample data as well as on synthetic data. We will describe how it is possible to derive from SR data not only the number of PA-FPs, as well as their photophysics, but also determine error bars around all of the estimates.

## Results

### Analysis of Simulated Data.

Consider the model shown in Fig. 1 for the photophysics of a single PA-FP. There are four possible states: inactive (*I*), active (*A*), dark (*D*), or photobleached (*B*). Once active, the PA-FP has two options: (*i*) it can blink possibly transitioning multiple times between the active and the dark state or (*ii*) it can irreversibly photobleach. Fluorescence is only detected from the active state, not in the other three states. Now consider a collection of *N* identical PA-FPs, each of which is governed by the model of Fig. 1.

Thresholding methods used in enumerating PA-FPs do not treat the stochastic nature of the fluorescence bursts explicitly. Thus, these methods cannot separate overlapping PA-FP signals where one FP activates before the last FP has photobleached. These methods also cannot treat superposing signals, which occur when two or more FPs fluoresce at the same time. Overlapping and superposing signals arise when multiple PA-FP spikes occur over a short period. This can happen when (*i*) there are a large number of PA-FPs; (*ii*) the PA-FP blinking rate is fast and the PA-FP photobleaching rate is slow; and (*iii*) the PA-FP photoactivation rate is fast and the PA-FP photobleaching rate is slow.

Without a priori knowing the PA-FP number, *N*, or the PA-FP’s properties, it is not possible to ascertain whether overlapping and superposing signals are improbable. Thus, in these conditions, if the problem at hand is specifically to determine *N*, the use of thresholding methods may not be appropriate. The method developed here, by contrast, is applicable for any *N* and PA-FP, regardless of its photophysical properties (activation, blinking, and photobleaching) because it explicitly treats them as the stochastic events they are.

To benchmark our method for finding kinetic rates and *N* from the data, we first analyzed 200 simulated trajectories, each one representing a collection of

These simulated trajectories (a sample of which is shown in Fig. S1) were generated using the Gillespie stochastic simulation algorithm (43). That is, the PA-FP photophysical dynamics were simulated stochastically according to the model given in Fig. 1 and the time trace recorded a fluorescence signal only during times in which any of the PA-FPs were in the *A* state.

There are several assumptions that these simulations make:

*i*) We assume that we have a unique dark and active state. If the number of dark and active states is greater than one, we will find that our estimates for the rates will come with large error bars (we will later see this for the case of Dendra2) and, in this case, additional information must be provided to specify the number of states needed.*ii*) We assume that the breadth of the distribution over*N*for our simulated data arises only from finite data. In reality, finite data provides a lower bound on the breadth of the distribution over*N.*The remainder of the distribution’s breadth arises from the intrinsic variability in the stoichiometry of the complex which we will discuss later.*iii*) We assume all PA-FPs photoactivate. Recent work suggests that the photoactivation efficiency ranges from only 40–60% depending on the PA-FP (44); to correct for this, we would amend our PA-FP topology to add another state to which state*I*can irreversibly transition without ever entering state*A.**iv*) We assume that maturation of all PA-FPs (mEos2 and Dendra2) used in the in vitro experiments is complete; we assume this because the elapsed time—from protein expression to the point where purified protein are obtained—is longer than the maturation half-time (90 and 120 min for Dendra2 and mEos2, respectively) (45, 46).*v*) We assume that protein complexes are well separated in space, which is often the case in experiments. For instance, the bacterial flagellar motor complex—with ∼22 MotB monomers forming part of each complex—are spatially well separated (20), as are individual kinetochores associated with six or more proteins linking kinetochore microtubules to centromeric DNA (47). For the complicated case where complexes are closely packed (i.e., packed within region spatially nonresolvable by PALM), our analysis should yield a distribution over*N*peaked at multiples of some number, e.g.,*N*= 5; this, in itself, is important information that would then motivate future experiments in an effort to establish whether complexes aggregate or coincidentally colocalize within an ROI.

We simulated three different blinking scenarios: moderate blinking (set 1), fast blinking (set 2), and slow blinking (set 3). We generated 200 synthetic traces for each scenario. The parameters for the three sets are summarized in Table 1. We tuned the ratio between the blinking rate and the photobleaching rate

To find the model parameters *Materials and Methods*. The likelihood function is based on the likelihood of observing a particular on-off fluorescence time trace in the data, which is then maximized with respect to the model parameters.

In the analysis of these data sets, we used a bootstrapping approach (resampling data with replacement) to determine the precision of our parameter estimates (48, 49). We randomly selected a subset of trajectories and determined the rates that maximized the sum of the log-likelihoods of the selected trajectories. We constructed a distribution of rates by repeating this process (200 times) for other randomly selected sets of trajectories. Our estimate for each parameter is the mean of the corresponding distribution, and we compute the 95% confidence interval based on percentiles of the distribution. The results are shown in Figs. 2 and 3 (for data set 1 and 2) and Fig. S2 (for data set 3).

The bootstrap results show that the parameters were determined precisely, except for a small underestimate of *N* throughout and a small bias toward slower *N* is due to missed transitions in the presence of fast blinking and low temporal resolution. In other words, as we increase the resolution to 5 ms, the bias in the estimate of *N* indeed vanishes.

Finally, Fig. 5 shows that if the values for the rates are known, we can use the PALM trajectory data to find the only unknown—namely, *N*. We show this to be the case in Fig. 5 for the challenging parameter set 1 both with 5-ms and 50-ms temporal resolution where the theoretical *N* value is set to 15. Again, the simulation with a resolution to 5 ms is sharply peaked at the correct answer.

### Analysis of in Vitro Data.

In addition to the simulated data presented above, we analyzed the in vitro data of Lee et al. (40). In this data set (which includes 1,000 traces), biotinylated Dendra2 molecules were immobilized on a streptavidin-coated glass coverslip, in a manner such that oligomerization was negligible and molecules were spatially well separated enough to be detected individually *I*) and then excite (from *I* to *A*) the Dendra2 molecules, respectively, until all of the molecules were photobleached (see Lee et al. (40) for further details). Individual emission bursts from the EMCCD output were processed into single-molecule time traces for analysis as described in ref. 40.

We simultaneously extracted the four kinetic rates (*N* distribution is sharply peaked at 1 (Fig. 6), as expected for this data set, because the experiments were designed to separate and isolate Dendra2 molecules on the coverslip. Our rate estimates (drawn from 300 bootstrap iterations) compare well with those of Lee et al. (40), who found a similar blinking rate (Table 2). Our results differ from their results most significantly for *D* to *A*. Lee et al. (40) fit the distribution of fluorescence-off times to determine

### Missed Transitions May Be Necessary in Tackling Larger Data Sets.

To tackle more challenging data sets, we need to correct for missed transitions when PA-FP residence times in the dark or active state are on the order of the temporal resolution or data acquisition time (the dead time, *Materials and Methods*. Typically, the lower bound for this time resolution is limited by the camera’s acquisition rate, i.e., the frame integration time.

In Fig. 7, *Left*, we show a distribution over *N* obtained by combining 10 traces *N* by a factor of ∼2 (even after 200 bootstrap iterations). This underestimation is most likely due to the missed transitions of the active state rather than of the dark state. Because the mean dwell time in the active state (83 ms

We did not undercount when we considered the

To approximately correct for the undercounting when *N* is larger than 1, we increased *N* estimate (Fig. 7).

Next, if we set the dead time closer to the active-state dwell time,

To avoid using this approximate treatment, a goal should be to develop PA-FPs that blink with kinetics slower than the data’s acquisition rate.

### Making Use of Intensity Measurements.

Previously, we showed how to correct for missed transitions in dealing with real data sets. Here we describe how to use fluorescence intensity data by separating bright states into several distinct aggregates also relevant to analyzing real data sets (Fig. 8).

Fig. 8*A* shows the number of states in each aggregated class for *Materials and Methods*), which determine the computational cost of computing the matrix exponentials and, ultimately, the computational cost of parameter extraction.

The method recovers *B*). A mathematical formalism is detailed in *Materials and Methods*.

Using the same number of synthetic traces and bootstrap iterations, we then extracted rates and *N* for a larger simulated data set *N* was underestimated because some transitions are missed. By regenerating simulated data with 5-ms time resolution, the rates are also correctly extracted and the estimate of *N* improves, but the distribution remains relatively broad because the data are probed to extract all rates as well as *N* (Fig. 9). The estimate for *N* sharpens dramatically, even when trying to estimate larger *N*, when information on the rates is provided and the data are instead entirely focused on estimating *N* more precisely (Fig. 5). All code used to generate the results can be obtained by emailing the corresponding author.

## Discussion

Common biochemical/biophysical methods used to establish the stoichiometry of protein complexes often lack single molecule resolution. For example, in the ion-channel literature the following methods are common: densitometry of silver-stained and Western-blotted protein complexes (in quantifying atrial inward rectifying potassium channel complex stoichiometry) (50); voltage-clamp electrical recordings [combined with channel inhibitors release to count potassium voltage gated channel (KCNE) subunits] (7); LC-MS (in counting subunit e used in bovine heart mitochondria ATP synthase self-association) (51); coimmunoprecipitation/immunocytochemistry (to assess colocalization of potassium voltage-gated channels such as KCNE2 to KCNQ1 and KCNE1 in cardiomyocytes) (17).

Fluorescence-based methods are also common in enumerating proteins but have severe limitations. For instance, flow cytometry (52, 53) yields estimates for intracellular protein number, but provides no spatial resolution and is less accurate for proteins numbering less than a few hundred (53). By contrast, quantitative fluorescence microscopy (54⇓–56), which is diffraction limited, provides a maximum spatial resolution of a few hundred nanomoles and is commonly used in protein counting, but still provides an order of magnitude less resolution than PALM or STORM (stochastic optical reconstruction microscopy), another SR technique. Using quantitative fluorescence, Joglekar et al. (54) counted proteins in the kinetochore of budding yeast using the fluorescence intensity from tagged Cse4p (57⇓–59) as a fluorescence intensity standard. This indirect method can lead to counting errors because standards themselves can be poorly characterized (18, 19).

Stepwise photobleaching methods, in turn, count proteins by monitoring discrete fluorescence intensity drops as individual proteins photobleach over one ROI. This method was used to count KCNE1s in the cardiac channel

By comparison, the method presented here provides a principled recipe for enumerating PA-FPs while gathering critical photophysical properties of the PA-FPs directly and self-consistently from the same SR data set used to determine *N*.

A qualitative reason explaining why likelihood functions peak at the correct value of *N* (Fig. S3) is illustrated with the following example. Suppose we start by guessing that all fluorescence spikes in a time trace are due to one PA-FP (i.e., *D* state, all waiting times between successive fluorescence spikes should therefore be sampled from this distribution: *N* value.

Another decisive advantage of our method is its ability to generate entire distributions of *N* for a given complex; this is particularly relevant given that protein complex stoichiometry may vary from complex to complex based on the dissociation constant of the various subunits in those complexes (5, 60), tissue development, and even its disease state (8, 61, 62). Developing therapeutic agents against some heart diseases, for instance, may require a quantitative assessment of protein complex composition in living tissues (7, 62⇓–64) and finding drugs to restore native complex stoichiometries is one such strategy (8, 62).

Our method is a first step in the quantitative determination of *N* from SR data, it includes the treatment of the “missed event problem” and shows how one can take advantage of intensity measurements.

## Materials and Methods

### States of the Model.

We assume that each fluorophore is independent of the others: the state of one fluorophore does not affect the state of any other fluorophore. We describe the state of a system of *N* fluorophores as a vector of the populations of the inactive, active, dark, and photobleached states: *I*, *A*, *D*, or *B*), and “macrostate” to refer to a population vector that describes the collection of fluorophores. For example, macrostate *i* for a collection of two PA-FPs in which both PA-FPs are inactive would be **S** is a population vector and *M* is the total number of macrostates.

Computing *M* is a common combinatorial problem: the number of unique ways to partition *N* indistinguishable objects into *x* bins. In this case, the objects are PA-FPs and the bins are the microstates. There are four microstates, so

For example, if the collection contains two fluorophores, the following 10 macrostates (obtained from

We model the collection of PA-FPs as an AMM. As shown in Table S1, each macrostate belongs to an aggregated class, either dark or bright. Macrostates with at least one active PA-FP

Here, the bright class corresponds to the detection of fluorescence and the dark class corresponds to the absence of fluorescence. The aggregated classes are necessary because PALM experiments cannot distinguish between the various dark states. The dark and bright classes here are analogous to the closed and open classes in the ion channel AMMs.

There can be multiple bright classes as there can be multiple levels of fluorescence that are observable: bright, bright2, bright3, …, brightN, as shown in Fig. 8.

The number of macrostates grows exponentially (Fig. S4) with *N*, which is a concern for the numerical calculations discussed in the results section; the computational time of the likelihood depends on the number of macrostates. As such, we define a quantity

### Model Kinetics.

The transition rate from macrostate

The dynamics of the PA-FP AMM are governed by a rate matrix, *t* are given by the Kolmogorov equation (34):

whose solution is given by

In the case of dark and bright observation classes, we can partition the rate matrix **Q** into four submatrices, based on the dark (subscript *d*) and bright (subscript *b*) aggregated classes:

The submatrix *d* to other states in class *d*; *d* to the states in class *b*. The other two submatrices are similarly defined.

Another possibility would be that multiple levels of fluorescence are observable: bright1, bright2, bright3, etc. For example, in the case of two bright aggregated classes, macrostates with

Having multiple bright aggregates is advantageous because it shrinks the size of the largest submatrix in the likelihood calculation. Furthermore, if the time traces infrequently visit the higher brightness classes, then fewer exponentiations of the largest submatrix will be necessary.

### Likelihood Function.

The likelihood function that we describe in this section provides an answer to the following question: Given a kinetic model and a set of kinetic rates, what’s the likelihood of observing the data? Our goal is to determine the kinetic rates and *N* that maximize the likelihood function, i.e., maximize the likelihood of observing the data that was collected by PALM.

Suppose we have a trajectory of *L* dwells representing the dynamics of a collection of *N* PA-FPs. Associated with each dwell is an observed aggregated class and a dwell time. The set of observation classes is *i* we observe class *t* and then transitioning to the bright class are given by the elements of the following matrix:

The *d* from its *d* for time *t*, and then transitioning to the *b*.

When we have missed transitions, we implement renormalized transition matrices—instead of transition matrices like those given by Eq. **5** —that account for missed transitions to the bright state, say, by resuming over all possible missed events (37, 65)

Next, we wish to calculate the likelihood of the dwell trajectory **h**, given the model parameters *θ*, where *θ* is the set of parameters (*N* and the transition rates) that go into the rate matrix **Q**. The likelihood function then reads as follows:

where *θ*, determine the elements of the rate matrix **Q**. The final factor of

### Numerical Evaluation of the Likelihood Function.

Our goal is to find the set of parameters **t** and **h**; in practice, this is accomplished by minimizing *N*, and then repeating this process for other values of *N*). The optimization is performed in Python using the limited-memory Broyden-Fletcher-Goldfarb-Shanno (BFGS) minimizer, implemented in Scipy (66). BFGS is a quasi-Newton method that performs well on nonsmooth optimization problems. The computational time is plotted in Fig. S5. The scaling depends on

As an aside, we point out that

To determine whether likelihood maximization runs converge to the correct parameter values, as a preliminary test we computed 1D slices in parameter space around the true value of each parameter. Each panel of Fig. S3 was obtained by holding four of the five parameters (*N*) constant at their true values and then varying the remaining parameter over a range encompassing its true value. We see that the four kinetic rates are peaked in the correct location at their true values. In the 1D slice for *N*, we see a maximum (as determined by the value of the likelihood function) that spans

Next, we assessed the ability of the numerical maximization procedure to converge to the correct parameter values. The 1D slices discussed above suggest that the likelihood function maximum is in the correct region of parameter space, but a separate question is, Can we simultaneously determine all five parameters? We found that the likelihood-maximization procedure converges to the correct kinetic rates within 100 cycles. Fig. S6 shows the convergence of the rate estimates from one of the maximization runs. For simulated fast (set 2) and slow blinking (set 3), we found that convergence of the likelihood maximization also occurred within ∼100 cycles.

## Acknowledgments

We thank S. H. Lee and A. Lee for providing the Dendra2 in vitro data and for numerous helpful discussions about PALM and the blinking problem. Support for this work was provided by a National Science Foundation (NSF) Graduate Research Fellowship (to G.C.R.); the Burroughs Wellcome Fund and NSF Grant MCB-1412259 (to S.P.); and the Howard Hughes Medical Institute and NIH Grant R01-GM032543 (to C.B.).

## Footnotes

- ↵
^{1}To whom correspondence should be addressed. Email: stevenpresse{at}gmail.com.

Author contributions: G.C.R. and S.P. designed research; G.C.R. and S.P. performed research; G.C.R., J.Y.S., C.B., and S.P. contributed new reagents/analytic tools; G.C.R. and S.P. analyzed data; and G.C.R., J.Y.S., C.B., and S.P. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

See Commentary on page 304.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1408071112/-/DCSupplemental.

## References

- ↵.
- Zamft B,
- Bintu L,
- Ishibashi T,
- Bustamante C

- ↵
- ↵
- ↵
- ↵.
- Nakajo K,
- Ulbrich MH,
- Kubo Y,
- Isacoff EY

- ↵.
- Osteen JD,
- Sampson KJ,
- Kass RS

- ↵.
- Morin TJ,
- Kobertz WR

- ↵.
- Yano M, et al.

- ↵.
- Lacoste TD, et al.

- ↵.
- Gordon MP,
- Ha T,
- Selvin PR

- ↵
- ↵.
- Qu X,
- Wu D,
- Mets L,
- Scherer NF

- ↵
- ↵.
- Wiedenmann J, et al.

- ↵
- ↵.
- Fron E, et al.

- ↵
- ↵.
- Coffman VC,
- Wu P,
- Parthun MR,
- Wu JQ

- ↵.
- Lawrimore J,
- Bloom KS,
- Salmon ED

- ↵
- ↵
- ↵.
- Hastie P, et al.

- ↵.
- Delalez NJ, et al.

- ↵
- ↵
- ↵
- ↵.
- Griesbeck O,
- Baird GS,
- Campbell RE,
- Zacharias DA,
- Tsien RY

- ↵
- ↵
- ↵.
- Lippincott-Schwartz J,
- Patterson GH

- ↵
- ↵
- ↵.
- Chudakov DM,
- Feofanov AV,
- Mudrik NN,
- Lukyanov S,
- Lukyanov KA

- ↵.
- Colquhoun D,
- Hawkes AG

- ↵.
- Kienker P

- ↵.
- Colquhoun D,
- Hawkes AG,
- Srodzinski K

- ↵
- ↵.
- Qin F,
- Auerbach A,
- Sachs F

- ↵
- ↵.
- Lee SH,
- Shin JY,
- Lee A,
- Bustamante C

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵.
- Johnston K, et al.

- ↵
- ↵.
- Efron B,
- Tibshirani R

- ↵.
- Corey S,
- Krapivinsky G,
- Krapivinsky L,
- Clapham DE

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵.
- Jiang M, et al.

- ↵.
- Bendahhou S, et al.

- ↵.
- Yano M, et al.

- ↵.
- Cingolani E, et al.

- ↵.
- Yu H, et al.

- ↵
- ↵