extfieldsofstudy:
Computer Science
Mathematics

text:
A/B Testing in Dense Large-Scale Networks: Design and Inference

Design of experiments and estimation of treatment effects in large-scale networks, in the presence of strong interference, is a challenging and important problem. Most existing methods' performance deteriorates as the density of the network increases. In this paper, we present a novel strategy for accurately estimating the causal effects of a class of treatments in a dense large-scale network. First, we design an approximate randomized controlled experiment, by solving an optimization problem to allocate treatments that mimic the competition effect. Then we apply an importance sampling adjustment to correct for the design bias in estimating treatment effects from experimental data. We provide theoretical guarantees, verify robustness in a simulation study, and validate the usefulness of our procedure in a real-world experiment.

Introduction
Measuring the effect of variants is a fundamental problem in several fields of study such as psychology (Rubin, 1974), epidemiology (Hernán & Robins, 2006), and many other sub-fields in medicine (Armitage et al., 2008). A/B testing is a commonly used method wherein randomized experiments are run with two or more variants. In the internet industry, A/B testing is very widely used in almost every aspect of product design and model improvements (Kohavi et al., 2013;2014;Tang et al., 2010;Xu et al., 2015).
Traditional A/B testing depends on the key assumption that the effect of treatment on an experiment unit is independent of the treatment allocation to other experiment units -commonly called "Stable Unit Treatment Value Assumption" (SUTVA) (Rubin, 1974). However, in many important network settings, this assumption is violated due to interference (Katzir et al., 2012;Toulis & Kao, 2013;Ugander et al., 2013). Estimating treatment effects under SUTVA violation is an important area of research. Several forms of network interference have been studied in the past (Backstrom & Kleinberg, 2011;Katzir et al., 2012;Aronow & Samii, 2012;Gui et al., 2015;Ugander et al., 2013).
In this paper, we focus on one such class of problems -a marketplace of commodity producers and consumers. Producers produce and make a commodity available for consumers, who in turn consume and give some desired utility back to the producers, which incentivizes them to keep producing. An example of such a marketplace is a social media platform, where the commodity is content and the return utility is feedback. Examples include Facebook, Instagram, LinkedIn, Twitter, etc. Other marketplace problems differ primarily on the commodity and return utility in question. Rides are commodities in Uber and Lyft, retail items are commodities in marketplaces like Amazon and eBay, and money is the common return utility.
In the marketplace setting, we define a treatment class represented by an edge-level boost factor, say T ij , where i is the producer and j is the consumer. From a consumer j's viewpoint, the boost factors T .j control or shape the consumer's exposure to different producers. This abstraction represents a very general class of treatments since any redistribution or shift in producer and consumer exposure can be expressed via appropriate edge-level boost factors. In content marketplaces, such treatments can reshape the content creator's exposure to their audience, whereas, in a transportation marketplace, it can reshape the likelihood of specific drivers and riders being matched.
For such treatments in marketplace problems, SUTVA is violated because of interference from the network. For instance, when we consider commodity producers as the experimental units, the effect of treatment on a producer is influenced by the effect of the treatment on all potential consumers of that producer, which in turn depend on each of their producers' experiences -thus leading to competition among producers connected to common consumers. For commodity consumers as well, SUTVA could be violated since a consumer's experience is determined by all of their potential producers.
When SUTVA is violated, then the primary principle used to measure the effect of treatment is to allocate experimental units as well as their first-degree neighbors (or first and second-degree neighbors, depending upon the interference arXiv:1901.10505v1 [stat.ME] 29 Jan 2019 function) to the same treatment (or control) (Eckles et al., 2017). The number of such units for which we can successfully allocate treatments in this fashion (i.e., whose neighborhood is appropriately treated) decrease with increasing density in the graph. Fewer experimental units result in lower statistical power in the measurements.
In this paper, we propose a novel technique named OASIS, "Optimal Allocation Strategy and Importance Sampling Adjustment", that provides a randomized testing framework for large-scale marketplace problems in the presence of interference. Our method actually works very well with dense graphs (with an increase in density being more helpful) and can be used to obtain high-power measurements as our approach does not require to allocate all first-degree (and second-degree) neighbors of an experimental unit to be in the same treatment. Our approach relies on the existence of an "intervening variable" -i.e., the effect of the treatment allocation to a unit's network on the experimental unit is fully captured by a sufficient statistic (cf. Arellano & Bover (1995); Cragg & Donald (1993);MacKinnon et al. (2002)). In marketplace problems, the total feedback received by a producer is a good example of such an intervening variable. The producers in a marketplace compete with each other for receiving feedback, and we attempt to mimic this competition effect in the design of experiment as much as possible. Furthermore, we construct an unbiased estimator of the treatment effect by applying an importance sampling correction. We show both from simulation studies and from a real-world experiment on a large social network graph that our technique works quite well for large-scale dense graphs. The method is also robust to small unknown deviations from the key assumptions. This technique can be used in any marketplace problem where the treatment in question can be expressed via the T ij abstraction, and where the key assumptions (outlined as Assumption 1, 2 and 3 in later sections) are satisfied.
The rest of the paper is organized as follows. In Section 2 we formulate the problem and discuss the setup in detail.
The key aspects of how we design an experiment via an optimization formulation are described in Section 3. We propose the OASIS estimator and provide theoretical guarantees in Section 4. In Section 5, we provide some empirical results both from an elaborate simulation study and a realworld experiment on a large social network graph, before concluding with a discussion in Section 6.

Problem Setup
A content marketplace can be represented as a graph, where every member is a node, and can be both a producer and a consumer. In mathematical terms, let us consider a network (or graph) G = (Ω, E), where Ω denotes the set of all members (or nodes) in G and E denotes the set of edges in G. We denote the set of neighbors of a member i in G by Ne(i) = {j : (i, j) ∈ E}. A producer i can only influence the set of consumers Ne(i) and a consumer j can only be influenced by the set of producers Ne(j). Let p base ij denote the normalized influence of producer i on consumer j so that we have i∈Ne(j) p base ij = 1 for all j. For example, p base ij can be the conditional probability that consumer j views a content produced by i given that consumer j views a content produced by Ne(j). A treatment T ij is a boost factor to the influence of a producer i on consumer j and we define T E = {T ij : (i, j) ∈ E}. Therefore, the normalized influence of producer i on consumer j in the treatment condition T E is given by (1) where Y (T E ) denote the response of a randomly chosen member when the entire population is under treatment condition T E . We consider m treatments T (1) E )} from a single randomized experiment. To this end, we randomly select m+1 disjoint subsets Ω 0 , . . . , Ω m of Ω. We aim to design an experiment T * E such that members in Ω r have the experience of the treatment condition T (r) E for all r = 0, . . . , m simultaneously. In order to quantify member experience, we make the following assumptions. Assumption 1. The response of member i depends on the treatment condition T E only through 1. The exposure of producers in Ne(i) to consumer i, defined as {p : k ∈ Ne(i)}, and 2. The total exposure of producer i to consumers, denoted by Z i (T E ).
Assumption 2. The total exposure of producer i to consumers Z i (T E ) depends on the treatment condition T E only through j∈Ne(i) α ij p for some constants α ij 's representing the strength of the relationship between producer i and consumer j. We further assume that α ij 's do not depend on the treatment condition T E .
An example of α ij can be the total number of interactions received by member i from member j when member j views i's content. In this example, Z i (T E ) would denote the total expected feedback received by member i. When α ij 's are known, it is natural to design an experiment T * E comprising of probabilities p * ij such that, 1. Exposure as consumer is satisfied. That is, for any i ∈ Ω r , r = 0, . . . , m and for all k ∈ Ne(i) we have 2. Exposure to consumer is matched. That is, If we are able to design such an experiment, then under Assumptions 1 and 2, an estimator of τ (T (r) In the following section, we discuss a strategy for designing T * E which satisfies the above two requirements under some constraints that control the risk and exposure of the experiment. We also theoretically demonstrate the robustness of our estimator under the violation of the assumption that α ij 's are known. In Section 4.2, we propose an importance sampling based modification of the estimatorτ (T (r) E ) that does not rely on the definition of the total exposure of a producer given in Assumption 2.

Design of Experiment
Let us denote Ω = ∪ m r=0 Ω r which we shall use for the rest of this paper. In this section, we construct an experiment T * E satisfying (2) and (3). Note that, since the probabilities are edge level, it is possible to easily obtain (2), by defining ki for all i ∈ Ω and k ∈ Ne(i). Thus, for all members chosen in the experiment, we can guarantee that their exposure as a consumer will match the true exposure as if the entire universe was exposed. Constructing T * E to satisfy (3) is a much harder task. We construct the experiment to satisfy (3) under the following constraints: 1. T * ij = 1 for all j ∈ S, where S is a randomly chosen subset of Ω satisfying E[|S|/|Ω|] = b and S ∩ Ω = ∅ for some pre-specified b < 1. This condition controls the risk of the experiment by forbidding a set of consumers S (disjoint from Ω ) to be exposed to treatment.

P[T
r∈0 Λ r is a set of consumers with predefined treatment allocation (described in details in Section 4.2). This condition controls the expected number of neighbors of a producer to be exposed to the experimental condition. The reason for letting q i 's depend on producer i is that we may want q i to be a function of the degree of member i.
3. T min ≤ T * ij ≤ T max , for some constants T min , T max satisfying 0 ≤ T min ≤ 1 ≤ T max . This third condition provides us a control over the maximum and the minimum value of the modified boost factors to be applied.
In Algorithm 1 we construct T * E , where for each producer i ∈ Ω we randomly choose a subset of neighbors ki for all k ∈ Ne(i), i ∈ Ω r ∪ Λ r and r ∈ {0, . . . , m}; 2: For each i ∈ Ω , construct a subset of neighbors C i by selecting each member in Ne(i) \ (Ω ∪ S ∪ m r=0 Λ r ) with probability q i ; 3: Obtain {p * ij : j ∈ C i , i ∈ Ω } by solving the optimization problem given in (5); C i respecting the first two conditions. Then we obtain {p * ij : j ∈ C i , i ∈ Ω } by solving a constrained optimization problem such that condition (3) is satisfied (detailed in Section 3.1). Finally, we obtain the modified boost factors as Note that that above choice of T * ij comes from the fact that if we use that in (1) we get, p

Optimization Problem
To satisfy (3) we would want to get p * ij which minimizes, under certain constraints. Recall that some of the p * ij are fixed by (2) and risk control (condition 1). Thus, we only need to find We add constraints to the optimization problem in a way such that T min ≤ T * ij ≤ T max . We choose 0 < R min < 1 < R max in order to control the ratio ∈ Ω } is defined as a solution of the following optimization problem:

Scaling the Problem
Solving the above optimization problem (5) is not an easy task especially when the number of edges n = |{(i, j) : i ∈ Ω , j ∈ C i }| is large. Even in moderately sized experiments in social network graphs, we can expect n to range in billions. Solving the generic quadratic programming (QP) problem in that scale is almost close to impossible. Thus, in order to be able to solve such a problem, we devise the following algorithm which solves it through an iterative approximation approach.
Let S denote the set of all consumers in the optimization formulation, that is, S = ∪ i∈Ω C i . We split this group of consumers into disjoint sets S k of roughly equal sizes, such that S = ∪ K k=1 S k . Since this induces a natural partition in the constraint space of (5), it is easy to see that any candidate solution {p old ij : i ∈ Ω , j ∈ S} of (5) can be improved by updating {p old ij : i ∈ Ω , j ∈ S k } with the solution of the following optimization problem, for all k ∈ {1, . . . , K}.
Note that (6) and (5) have the same set of constraints for each j ∈ S k , since {i : We provide an efficient strategy in Algorithm 2 to solve the overall optimization problem. We start with a feasible candidate p old ij = p base ij and run an iterative scheme to update p ij using (6) as we loop over each k = {1, . . . , K}. Once this inner loop completes, we get the full next best {p ij : i ∈ Ω , j ∈ S}. We continue the outer loop till convergence. By doing this iterative scheme we are able to solve much larger problems, since size of the each optimization problem n k := |{(i, j) : i ∈ Ne(j) ∩ Ω , j ∈ S k }| is much smaller than n. In fact, the worst-case complexity of the iterative method with nIter outer iterations is given by O(nIter × K k=1 n 3 k ), while the worst-case complexity of the original optimization problem (5) is O(( K k=1 n k ) 3 ). Moreover, we can tune on the size of the partition of S such that we can optimize the on the total time given the memory restriction on the QP solver. We use the Operator Splitting method to solve individual QP (Stellato et al., 2017;Banjac et al., 2017a;. Algorithm 2 Solving for Optimal p ij 1: Input: Group Size g, tolerance level tol, Outer Iteration Limit maxIter 2: Initialize the p ij = p Solve the optimization problem as in (6) ij for all (i, j). 10: end for 11: Return the optimal p * ij = p (maxIter+1) ij

Theoretical Analysis
In this section, we theoretically demonstrate the robustness of our estimator under the violation that the α ij 's are known. We also derive an unbiased estimator based of τ (T (r) E ) via importance sampling that does not depend on the definition of total exposure. All proofs are pushed to the supplementary materials for brevity.

Robustness
In order to demonstrate robustness of our experimental design when Assumptions 1 and 2 hold but α ij 's are not known, let us assume that we construct an experiment T * E by randomly perturbing p for some known β ij 's. In the following theorem, we show that as long as the random perturbation does not depend on β ij p where g is a differentiable function with respect to the first coordinate, such that {U * ij : j ∈ Ne(i)} are identically distributed random variables satisfying If there exist an experimental design T * E such that

Importance Sampling Adjustment
We have assumed that response of producer i depends on the total exposure or feedback obtained from all consumers which we denoted by Z i (T E ). That in turn was assumed to depend only through j∈Ne(i) α ij p . As a result, we wanted to match on the above quantity. However, in many cases, we would not be able to exactly match this quantity for every i. As a result, we only requested an approximate matching (3) in our experimental design. In this section, we describe a method that allows us to get an unbiased estimator of τ (T (r) E ), even when we are only able to get an approximate matching. We begin with a formal assumption. Assumption 3. There exists an observable random variable X i (T E ) such that the conditional distribution of Y i (T E ) given X i (T E ) depends on the treatment condition only through W i (T E ), an unknown function of {p Note that when Assumption 1 holds, Assumption 3 implies that Y i (T E ) and Z i (T E ) are independent conditionally on X i (T E ) for any treatment condition T E . Furthermore, Theorem 2. Let f r denote the density of X(T (r) E ), and let f * r denote the density of X i (T * E ) conditionally on i ∈ Ω r , where T * E is the output of Algorithm 1. We definê .
Remark. Note that this allows us to skip Assumption 2 while doing the experimental design. However, from a practical standpoint, if we (approximately) know the dependency of Z i (T E ), we should definitely use it to do the matching while designing the experiment. The design of experiment step is crucial for reducing the variance of the importance sampling weights as well as the corresponding estimator.

Density Estimation
The densities f * r and f r are unknown and we need to estimate them from the data. Since we observe {X i (T * E ) : i ∈ Ω r }, any parametric or non-parametric density estimation method can be applied for estimating f * r . For estimating f r , we make the following additional assumptions: , where X 1 d = X 2 denotes that X 1 and X 2 have identical probability distributions. E ) : j ∈ Ne(i) ∩ Λ r , i ∈ Ω r , r = 0, . . . , m} are observable for sets of members Λ r ⊆ ∪ i∈Ωr Ne(i) (see Algorithm 3).

Assumption 5. For each treatment condition T
Assumption 4 allows us to obtain an estimate of f r by adjusting the mean and the variance of an estimated f * r . Assumption 5 facilitates the estimation of moments of X(T (r) E ) from the data {X ij (T (r) E ) : j ∈ Ne(i) ∩ Λ, i ∈ Ω r }. In following theorem, we propose consistent estimators of the first and the second moments of X(T (r) E ). Theorem 3. Fix r ∈ {0, . . . , m}. Let Λ r be as in Assumption 5. Assume that are independent of |Ne(i)|, µ i and σ 2 i . We define Then, under Assumption 5, In the following algorithm, we propose a method for selecting {Λ 0 , . . . , Λ m } such that Λ r ⊆ ∪ i∈Ωr Ne(i). It is easy to verify that Algorithm 3 does not introduce any selection bias for choosing Λ r in the sense that each member in ∪ i∈Ωr Ne(i) has the same probability to be included in Λ r . Note that we use the same {Λ 0 , . . . , Λ m } in Algorithm 1.

OASIS
We conclude this section by gluing all the pieces together in the following algorithm, called Optimal Allocation Strategy and Importance Sampling Adjustment (OASIS).

Algorithm 4 OASIS
Output:τ (T Computeτ (T (r) where Φ is the CDF of the standard normal distribution.

Experiments
In this section, we describe an in-depth simulation study to demonstrate the robustness of our method under moderate violations of Assumptions 2 and 3. Furthermore, we compare OASIS with an oracle cluster-based method that uses a given set of clusters for applying full treatments to all members in a randomly chosen cluster. We also apply our technique to a real-world experiment on the LinkedIn graph. We were able to validate our mechanism by being able to match the results for a full population experiment. The details are given below.

Simulation Study
We generate a random graph G with |Ω| = 50000 vertices and the average degree equals 100 as follows. We first obtain a graph G BA = (Ω, E BA ) by combining 10 randomly generated graphs with 5000 vertices and average degree equals 80, where each graph is generated according to the Barabasi-Albert model (Barabási & Albert, 1999) with the power of the preferential attachment equals 0.25. Then we generate an Erdös-Rényi graph (Erdös & Rényi, 1959) G ER = (Ω, E ER ) with 50000 vertices and average degree equals 20, where all pairs of nodes have an equal probability of being connected. Finally, we obtain G = (Ω, E BA ∪ E ER ). For each ordered pair of edge (i, j), we independently generate α ij and p base ij from the following distributions: where d j denotes the degree of node j, U ij 's are i.i.d. U nif orm[10, 100] random variables and V ij 's are i.i.d. U nif orm[1, 2] random variables. For all treatment conditions T E , we define where η ij 's are independently distributed N (0, (d i d j ) −1 ) random variables. Next, we define where g(x) = 10/(1 + exp(−x/10)), β ∈ {0, 1} and i 's are independently distributed as a standard Gaussian random variable. Note that β = 0 corresponds to the absence of the second order effect. Furthermore, the strength of the second order effect Z i (T E , δ) decreases as δ increases. Finally, we define We consider three different treatment conditions corresponding to γ ∈ {0, 0.5, 1.5}, where γ = 0 corresponds to the control case.
It is easy to see that our simulation setting respects Assumption 1. Assumption 2 is satisfied if and only if δ = 1 and α ij 's are known, and Assumption 3 is satisfied if Z ij (T E , δ)'s are observable. In order to verify robustness of our method, we do not assume that α ij 's are known, and we apply Algorithm 1 for constructing T * E assuming α ij = 1. Furthermore, we assume that we are able to observe only a noisy version of Z ij (T E , δ), given by Therefore, we apply the importance sampling correction based on estimated densities of X i (T * E , δ) and X i (T E (γ, δ)), where X i (T E , δ) := j∈Ne(i) X ij (T E , δ). In particular, we assume that both X i (T * E , δ) and X i (T E (γ, δ)) have Gaussian distributions and we estimate the moments using the technique described in Section 4.3. The detailed parameter setting are given in the supplementary material.
We use OASIS (see Algorithm 4) with the self-normalized importance sampling 1 wi Y i w i , as it is known to have a more stable behavior in practice. Furthermore, we apply the bootstrap method for estimating the variances of the OASIS estimators, as described in Section 4.4. Table 1 demonstrates that the OASIS estimator can achieve satisfactory coverage probability even under violations of the underlying assumptions, except for the few cases corresponding to strong violations of assumptions 1 .
Next, we compare OASIS with an oracle cluster-based method that takes advantage of the presence of disjoint connected components in G BA to design the following experiment: where H(γ)'s are randomly chosen disjoint connected components of G BA for γ ∈ {0, 0.5, 1.5}.   Figure 1 shows that OASIS performs equally well or slightly better than the oracle cluster-based method in all simulation settings, except for β = 1 and δ = 1.5. This could be explained by the facts that (i) smaller values of δ tend to induce more bias in the cluster-based estimator, and (ii) larger values of δ lead to stronger violations of Assumption 3 as the observed variable X i (T E , δ) suffer from low signalto-noise ratio.

Real World Experiments
Although the above-described simulation study shows the robustness of our technique, we were interested in validating the result on a real-world experiment. In order to do so, we considered the LinkedIn social network graph. For simplicity, we only considered a single experimental setting T E = {T ij : (i, j) ∈ E}, where T ij are the boost factors for boosting content from producer i to consumer j and E is the set of all edges in the LinkedIn connection graph. We aim to validate our method by matching the feed metric measurements with a proxy of the 100% ramp on Ω. Since we lack an oracle, in this case, we create the following group of randomized members.
1. True Treatment Producer (M 1 )-These members have their entire first degree exposed to the treatment. Moreover, the first degree will see boosts from all eligible producers. Their behavior will mimic the 100% ramp on Ω. We fix this to 0.5% of the members.
2. Balanced Treatment Producer (Ω 1 ) For these members, we apply our mechanism. That is, we randomly pick some of their first degree members and appropriately change the boost factors so that the total interaction received (an approximation of Z i (T E )) is similar in distribution to M 1 . We keep this set to 1%.
This experimental setup is slightly different than that usual A/B testing setup. Here, we want to match the OASIS estimate with an unbiased estimate of the true effect, obtained from M 1 . Note that we cannot have a pure control group simultaneously with M 1 , since the graph is not separable. Thus, we compare the treatment effects with the control effect estimated before the start of the experiment. We would consider this experiment a success if we are able to match the metric reads from M 1 with Ω 1 . Remark. Even with such a small ramp, the total number of members and edges that were exposed to T * ij = 1 was about 13 million and 130 million respectively.
We ran the experiment for a week and observed that there was no significant impact on the consumer side metrics (such as clicks, likes, comments, etc). For the producer side metrics (such post creations, contributions, etc) we saw that there were significant differences between both M 1 vs control and Ω 1 vs control, which was expected. For proper validation, as explained before, we were interested in comparing M 1 vs Ω 1 . In doing so, four metrics came to be significant at the 5% level, namely mobile messaged from feed (m 1 ), mobile shares from feed (m 2 ), all messages from feed uniques (m 3 ) and mobile messages from feed uniques (m 4 ). However, after importance sampling adjustment, we were able to correct the p-values and all metrics were insignificant at the 5% level (see Table 2). All other feed related producer side metrics when compared between M 1 and Ω 1 were non-significant, thus verifying our network measurement mechanism was working as expected in real dense graphs. All parameter settings are given in the supplementary material.

Discussion
We have presented a two-step method, called OASIS, for estimating the average treatment effect for a class of continuous treatments in networks with interference. First, we design an experiment by optimally allocating treatment exposure to a set of randomly selected consumer-producer pairs in the network. We solve a large-scale linearly constrained quadratic program to achieve this optimal allocation of the treatment exposure. Secondly, we apply an importance sampling correction for estimating the average treatment effect from the experimental data. The importance sampling adjustment corrects for the design bias induced by the violations of assumptions and/or the restrictions applied to the optimization problem for controlling the risk of the experiment. We estimate the variance of the OASIS estimator by applying the bootstrap method. The correctness of the OASIS estimator relies on a number of assumptions, and we demonstrate the robustness of OASIS to moderate violations of these assumptions with an extensive simulation study. Finally, we validated the usefulness of OASIS with a real-world experiment on a large network.
OASIS provides a number of interpretable, tuning knobs for controlling the risk of the experiment as compared to cluster-based methods. These knobs can be set to their extreme values to make the design the most favorable to the OASIS estimator. A crucial advantage of the OASIS estimator is that it tends to perform better for dense networks, while the cluster-based method would have an advantage over OASIS for sparse networks that can be easily decomposed into clusters. An interesting future work could be to combine a cluster-based approach with OASIS in order to gain additional robustness and efficiency.
For some marketplaces (e.g., Uber, Lyft, Amazon), the graph is inherently dynamic or partially known. Our work assumes that the graph is completely known and fixed. Approximating a dynamic graph with a static graph which includes the various possibilities can work, but would not scale well. Hence, we will need to extend our method to partially known and/or temporally dynamic networks. We recognize this as an important and challenging future work.
This completes the proof. Therefore, where the last equality follows from (9) and the fact that ρ i and µ i d i are independently distributed.
Next, it follows from similar calculations that Furthermore, it is easy to verify that ρ 2 Using the independence of (ρ i , ρ i ) and (d i , µ i , σ i ), we obtain Hence, we have,


Recovering the Imperfect: Cell Segmentation in the Presence of Dynamically Localized Proteins

Deploying off-the-shelf segmentation networks on biomedical data has become common practice, yet if structures of interest in an image sequence are visible only temporarily, existing frame-by-frame methods fail. In this paper, we provide a solution to segmentation of imperfect data through time based on temporal propagation and uncertainty estimation. We integrate uncertainty estimation into Mask R-CNN network and propagate motion-corrected segmentation masks from frames with low uncertainty to those frames with high uncertainty to handle temporary loss of signal for segmentation. We demonstrate the value of this approach over frame-by-frame segmentation and regular temporal propagation on data from human embryonic kidney (HEK293T) cells transiently transfected with a fluorescent protein that moves in and out of the nucleus over time. The method presented here will empower microscopic experiments aimed at understanding molecular and cellular function.

Introduction
In the past decades, it has become evident that proteins are very dynamic and their localization within the cell dictates which function they perform [3]. Cell biologists carry out time-lapse fluorescence microscopy experiments to study protein localization and unravel how this affects the protein's function. Optogenetics can be used to control the localization of the protein of interest [18] to observe the effects of dynamic localization patterns. Take Figure 1 for an example: HEK293T cells expressing mCherry fused to the optogenetic tool LINuS [15] were observed over time to analyze the effect of blue light on the nuclear import of the protein. In the absence of blue light, the fusion protein localizes predominantly in the cytosol. Hence, the nucleus appears dark. However, in presence of blue light, the fusion protein enters the nucleus making it appear bright and difficult to distinguish from the cytosol. Giving light repeatedly creates oscillations of the protein in and out of the nucleus in time. Learning-based methods, such as U-Net [19] and Mask R-CNN [6], succeed at segmenting structures in data with clearly visible patterns, but fail when the visibility deteriorates. When the signal in the nucleus is similar to that in the cytosol, nuclei segmentation from any network is not reliable. While in a single image it is not possible to improve these segmentations, past and future frames in a video provide additional information for refinement. An expert can play the video back and forth to infer the segmentation of ambiguous nuclei. An automated segmentation method, too, must (1) automatically identify critical frames and (2) propagate predictions from neighbouring frames.
In this paper, we address both challenges by (1) equipping Mask R-CNN with uncertainty estimation to identify erroneous predictions and (2) incorporating optical flow to improve the identified erroneous predictions by propagating certain predictions from neighboring frames. Doing so, we introduce the most recent uncertainty estimation methods in biomedical instance segmentation and solve a real task commonly experienced in signalling studies which is not yet addressed. So far, nuclear markers have been employed in the experiments so that the available automated segmentation tools can be used [4]. However, additional markers cause unreliable quantification since different proteins bleed-through and interfere with each other. They also limit the channel space needed for other proteins of interest. The presented method makes the use of a nuclear marker dispensable.

Related Work
We are the first to address instance segmentation of structures with an oscillatory fluorescent signal in biomedical videos. Our work is related to methods that benefit from temporal features in their design. Milan et al. [14] and Payer et al. [17] used Recurrent Neural Network (RNN) to aggregate temporal features. Paul et al. [16] incorporated temporal cues for segmentation via optical flow. Similarly, Jain et al. [9] propagated features instead of segmentation with flow from key frames. A similar idea was applied to instance segmentation by Bertasius and Torresani [2]. Although these methods seem close, the task is different: (1) they heavily rely on dense annotation in time to learn interpolations explicitly while we cannot afford it due to high cost of expert annotation and difficulty in fine-grained annotation of imperfect frames and (2) they benchmark only on visible objects, while we are solely interested in objects with limited visibility.

Instance-Aware Segmentation with Uncertainty Estimation
We base our model on Mask R-CNN [6] equipped with elastic deformations of U-Net [19] to create additional biomedically plausible images on-the-fly for better generalization. We incorporate uncertainty estimation in our Mask R-CNN architecture to detect erroneous predictions; see Figure 2. We consider data uncertainty (aleatoric), model uncertainty (epistemic), and their combination.
Data uncertainty. For data uncertainty, we use the modified cross-entropy loss [10] in the mask branch of Mask R-CNN. This models the data uncertainty as the learned noise scale from the data. To learn both the class scores and their noise scale, the negative expected log-likelihood is minimized for pixel i as: where c is the correct class among all classes (c) andŝ are predicted logits corrupted by Gaussian noise with standard deviation σ. σ is also learned by the network alongside with the logits. Figure 2 shows where this loss is embedded in the network. At test time, the uncertainty is computed as the entropy of the class pseudo-probabilities for each pixel as: Model uncertainty. Model uncertainty requires sampling from the model. We experimented with the latest sampling strategies. Dropout is one of the common techniques to sample from networks [10]. Ensemble replaces the sampling

Algorithm 1: Uncertainty-Based Nuclei Propagation
Input: instance segmentation masks S, uncertainty maps U, motion per consecutive frame pairs M, hard threshold θ, relative thresholds: α,β Output: tracks for frame id f in order do by dropout by individually trained networks [11]. SGDR ensemble replaces the individual trained networks by the pre-converged snapshots of the same model [7]. Pre-converged models are obtained at the end of each SGDR (SGD with warm restarts) [12] cycle. WTA (Winner-Takes-All) [8] aims for training a single network with multiple heads where at each iteration only the head with the best prediction gets penalized. EWTA (Evolving WTA) [13] is a variant of WTA which improves the trade-off between diversity and consistency.

Uncertainty-Based Nuclei Propagation
We base the nuclei propagation on the cell tracker 4 by Ronneberger et al. [19]. This overlapping-based tracker is well suited for the data at hand, which does not show large motion over time. Despite its simplicity, it is one of the top performing cell trackers in the challenge. After computing the overlapping-based tracks, we go over all the frames of each track with ascending average nuclei uncertainty order (without nuclei, uncertainty is set to infinity) and update the identified uncertain predictions by propagating masks from more certain neighbours as depicted in Figure 1. A prediction is marked uncertain if its average uncertainty is higher than a threshold θ. Neighbours are more certain if they meet a relative threshold of α in case of a single-side propagation and β in case of a two-sided propagation as detailed in Algorithm 1. This resembles the procedure that experts follow to find non-visible nuclei by sliding over time.

Motion Estimation for Biomedical Videos
Propagation of predictions over time requires motion estimation between frames to warp the certain predictions onto the less certain ones. One simple way is to warp by the shift and scaling parameters computed between the not yet updated and neighbouring nuclei predictions. This approach assumes the shape of the nuclei does not change over time; however, slight deformations can occur. Optical flow can provide fine-grained motion. Recent optical flow methods for natural images are networks [5] trained on synthetic datasets. These methods perform well on real images, but their performance deteriorates as the gap between real and synthetic images grows. Our data is very different from existing synthetic datasets and no synthetic data exists for biomedical data. Therefore, with the prior knowledge that expected flow in our videos can be well explained by smooth deformations, and being able to generate them on the fly via deformation augmentations of U-Net, we explicitly train a network to predict these deformations as explored by Sokooti et al. [20]. We warp an image (I t+1 ) backward in time with randomly generated smooth deformations (f t→t+1 ) to obtain the previous image (I t ). Then we train a FlowNet on the on-the-fly-generated image pairs and ground-truth flows. Even though the fluorescent signal in the nuclei regularly disappears making the optical flow challenging to estimate, the motion of the cytosol helps infer the motion of the nuclei.

Experiments
Implementation details. We based our implementation on the Mask R-CNN by Abdulla [1]. We incorporated publicly available elastic deformations 5 for U-Net [19] as additional augmentation. We used ResNet50-FPN as backbone. We used softmax cross entropy for the mask head with 3 classes (background, cytosol and nuclei). We trained the networks from scratch for 12k epochs except the SGDR ensemble, which was trained for 15k with 3k cycles. We used the last pre-converged models. The WTA and EWTA models were trained for 11k and then merged with a network trained for 4k. We used 4 as ensemble size for all methods. Ensembling was performed over bounding boxes. For nuclei propagation, we used θ = 0.5, α = 0.7 and β = 0.85. For Mask R-CNN and FlowNet (variant C) training we used deformations with 3 and 10 control points, respectively and deformation magnitude of 10. Please see the original papers for more details. Our code and examplar dataset is publicly available 6 . Data and annotation. While the challenge we address is common in signalling studies, there is no public dataset for this purpose. Our data was generated after 24 hours of transient transfection of human embryonic kidney (HEK293T) cells with a construct expressing the fusion protein mCherry-LINuS with a ZEISS LSM780 confocal microscope. MaskRCNN was trained on 82 cropped images from original 1800x1800 images with pixel size 0.11x0.11 and 6 cells on average in each and tested on 2 full-sized images with 73 cells in total to validate instance segmentation and uncertainty estimation. We evaluated our propagation algorithm on an unseen video with 35 full-sized frames with randomly selected 117 cells annotated by experts. In cases where the nuclei appeared ambiguous experts interpolated the annotations from past and future frames.

Instance Segmentation and Uncertainty Estimation Evaluation
We evaluated our instance segmentation by mean average precision (mAP). In Table 1 we refer to standard mAP as mAP (sm) since it is based on softmax scores. To further evaluate the quality of the uncertainty estimation, we replaced the softmax scores by average entropy over the cell and re-computed mAP (mAP (ent)). This simulates our approach as we rely on averaged uncertainties in our nuclei propagation. This measure shows how reliable the predicted uncertainties are at ranking the prediction quality in the precision-recall curves, which are used commonly to evaluate uncertainty estimation for classification tasks. In Table 1 we see that WTA merged is the best at 0.5 IoU (Intersection over Union) threshold and at 0.75 threshold ensemble is the best. In the rest of our experiments we used the WTA merged with data uncertainty since it is computationally more efficient than ensemble.

Nuclei Propagation Evaluation
We report the mean IoU of the nuclei segmentation for all experiments in Table 2. We use interpolated for nuclei completely missed before our improvement, updated for the nuclei that were segmented but our method decided to improve, non-updated for the nuclei that were chosen not to be improved by our method and all for all the mentioned cases. The first row shows the results before any propagation algorithm. We present 3 variants of motion used in propagation (column: warped with). We also explored using the certainties as pixel-wise weights in fusing segmentations from candidate neighbors and computed the weighted average to find the final mask. To isolate the gain by our uncertainty-based error detection, we created a baseline which is identical to our mean-flow variant, but performs the propagation on all the nuclei over a track. The significant improvement obtained by all our variants, especially flow variants without fusion, shows that our method can effectively improve erroneous nuclei predictions. The baseline that propagates to all frames (all ) had a lower performance on the non-updated frames, showing that propagation independent of an uncertainty measure harms nuclei with high confidence (approx. 90% of all nuclei).

Slight Signal Loss
Since the network has been trained on a very clean training data, even small shifts in the signal distribution might disturb its performance and cause failures as in Figure 3. Notice that in the first frame, we have two connected components wrongly building together the nucleus mask. As each cell can only have a single nucleus, we consider propagating one which is the most certain. The remaining connected component can be filtered with a simple post-processing step. We see that our method improves upon the initial segmentation masks significantly.

Extreme Signal Loss
The scenario in Figure 4 shows the full potential of our method in the case of extreme signal loss. Uncertainty estimation allows to rank segmentation masks and infer more reasonable masks for the less certain frames. In the last frame the network misses the nucleus completely in a very confident way. Since we set the uncertainty in such cases to infinity to ensure that it gets updated (biologically there is certainly a nucleus in each cell), we are able to recover the missing nucleus. For a biologist taking measurements from the results of our method is more reliable compared to the non-updated results. Fig. 3: Example of a nucleus getting slightly less visible due to increased noise at time t. Columns correspond to time points t − 1, t and t + 1, from left to right respectively. In the first row we show the raw images. In the second row we overlay MaskRCNN outputs with the input image and report the average uncertainty. Pixel-wise uncertainty maps are presented in the third row while updated masks are overlaid to the original images in the last row.

No Signal Loss -Segmentation Failure
If the training data is not rich with all variations of the real task at hand, which is often the case in biomedical settings, networks can also fail at generalization. Our method is not only limited to refining failures due to signal loss, on the contrary, it is generic to all failures caused by the uncertainties inherited in the network or the data. Our predicted uncertainty estimations provide a good indicator for these scenarios and allows to recover them as we show in Figure 5.

Conclusion
We addressed automated segmentation of image sequences, which cannot be analyzed frame-by-frame due to temporary uncertainties causing errors in predictions. First, we estimate uncertainty from Mask R-CNN to identify unreliable predictions. Second, we improve the less reliable predictions by propagating Fig. 4: Example of a nucleus becoming completely indistinguishable from the surrounding cytosol in all time points. Columns correspond to time points t − 1, t, t + 1 and t + 2, from left to right respectively. In the first row we show the raw images. In the second row we overlay MaskRCNN outputs with the input image and report the average uncertainty. Pixel-wise uncertainty maps are presented in the third row while updated masks are overlaid to the original images in the last row.
the more certain ones from neighbouring frames. We evaluated our method on HEK293T cells expressing a protein that oscillates in and out of the nucleus over time making the nucleus invisible temporarily. Our method improves nuclei segmentation over several baselines while keeping the cytosol and background segmentation untouched. We believe that our method will facilitate further timeseries analysis for quantitative biology to understand the effect that the dynamic localization of the protein has on the cell without additional markers. Columns correspond to time points t − 1, t and t + 1, from left to right respectively. In the first row we show the raw images. In the second row we overlay MaskRCNN outputs with the input image and report the average uncertainty. Pixel-wise uncertainty maps are presented in the third row while updated masks are overlaid to the original images in the last row.


extfieldsofstudy: Physics Medicine

text:
A density-based threshold model for evaluating the separation of particles in heterogeneous mixtures with curvilinear microfluidic channels

Particle separation techniques play an important role in biomedical research. Inertial focusing based microfluidics using nonlinear channels is one of the promising label-free technologies for biological applications. The particle separation is achieved as a result of the combination of inertial lift force (FL) and Dean drag force (FD). Although the mathematical expressions of FL and FD have been well derived in prior studies, they are still complicated, which limits their popularity in practice. Recent studies modified these expressions through experiments and proposed a threshold model, which assumes that only particles larger than the threshold will be well focused. Although this threshold model has been used in recent studies, two varying versions of the threshold model (TM1 and TM2) prevents standardisation in practice. In addition, both models were developed with regular low-density particles and may not be applicable to samples with higher density or samples with irregular shapes. Here, we evaluated the threshold models with samples of different densities. Based on these evaluations, we derived a modified model (TM4), which additionally considers the factor of particle density to improve the accuracy of existing models. Our results demonstrated that TM4 could more reliably predict the sorting efficiency of samples within a wider density range.

Scientific Reports | (2020) 10:18984 | https://doi.org/10.1038/s41598-020-75878-w www.nature.com/scientificreports/ where ρ is the fluid density, U is the maximum fluid velocity, μ is the fluid viscosity, and D is the hydraulic diameter expressed as: where h and w are the height and width (for rectangular cross-section), respectively. In inertial microfluidics, the resultant force due to shear-gradient and wall-induced forces acting on particles in the opposite direction is called F L, which is given by: where C L is lift co-efficient, which is a function of particle position, and a is the particle diameter. Apart from inertial lift force, particles flowing in the curved channel will experience a F D resulting from 2 counter-rotating vortices (Dean vortices), which can be quantified by Dean number (De) as: where R is the radius of curvature of the channel path. The Dean vortices induced at the channel cross-section will cause particles to move back and forth across channel width. This lateral velocity for particles, Dean velocity (U D ), is estimated by: With Dean velocity (U D ), the F D acting on the particles can be expressed as: The focusing effect for given particle dimensions is determined by the ratio of F L and F D (R f = F L /F D ), which varies with particle diameter (R f ∝ a 3 ). The largest particles will be focused on the channel position closest to the inner wall as they experience the larger F L compared to other particles at that position in the system. A threshold model was derived based on the modification of previously established Eqs. (3), (4), (5), and (6), among which only particles larger than the threshold will experience significant F L, thereby balancing F D to facilitate particle focusing. The threshold of particle focusing is proportional to the size of the microchannel, which can satisfy the following inequality, as stated: where x = D( hydraulic diameter) for TM1 and x = h (channel height) for TM2.
Particle recovery rate (PRR) is often used to quantify the particle sorting performance in microfluidics 23 . In this study, the PRR of the device used ( Fig. 1) is defined as: Evaluation of TM1 and TM2 with regular low-density particles. Given that TM1 and TM2 were developed with regular low-density particles ( ∼ 1.05 g/ml), both threshold models were initially analysed by COMSOL, and then polystyrene (PS) microsphere experiments were performed. According to Eq. (7), the estimated thresholds of TM1 and TM2 are 21 μm and 15 μm, respectively. Therefore, 21 μm and 15 μm regular PS particles were first studied. In the COMSOL simulation, it was found that 21 μm particles (threshold of TM1) were collected at outlet 4 ( Fig. 2a), and 15 μm particles (threshold of TM2) were collected at outlets 3 and 4, respectively (Fig. 2b). Using the estimated values provided by COMSOL, we proceeded to conduct actual experiments to validate this hypothesis further.
For the experiments with 15 μm and 21 μm PS microspheres, more than 90% (91.6 ± 7.3%) of 21 μm particles were recovered in outlet 4 (Fig. 2c), and the rest of the particles was found at outlets 3 and 5 (5.2 ± 3% and 3.2 ± 4.4%, respectively). The particles found in outlets 3 and 5 could be a result of unfocused streamlines in the first few seconds of flow (Fig. S1). For 15 μm particles, the PRRs at outlets 3 and 4 were 52.4 ± 4.6% and 47.6 ± 4.7%, respectively (Fig. 2d). The highest PRR of 21 μm microspheres indicated that the threshold estimated by TM1 for PS microspheres was more accurate than TM2 for devices with similar spiral radius and aspect ratio parameters.
Evaluating the overestimation of TM1 threshold. The previous results cannot exclude the possibility of overestimating the threshold. For example, even if the threshold in the system was 21 μm (TM1), particles larger than 21 μm (such as 25 μm) would still be well focused, but 25 μm was definitely not the precise threshold value. To address this issue, we studied the PRRs of microspheres in the range of 21 to 15 μm. If the recovery     (Table S2). All 18 μm particles (T t1 = 18 μm) were well focused on outlet 4 ( Fig. 3a), which was similar to the PRR of 21 μm particles (Fig. 2a). Hence, the threshold value of 16.5 μm (T t2 = (18 + 15)/2 = 16.5 μm) (test 2) was used for the next test. Then, all 16.5 μm particles were well focused on outlet 4 ( Fig. 3b), which was similar to the PRR of 18 μm particles (T t1 ). Hence the threshold value of 15.7 μm (T t3 = (16.5 + 15)/2 ∼ 15.7 μm) was used for the next test (test 3). The simulation of 15.7 μm particles showed that only 66.7% particles (n = 3) were concentrated on outlet 4 (Fig. 3c), significantly lower than the PRR of 16.5 μm particles. As a result, the target threshold range was likely to be between 16.5 μm and 15.7 μm (15.7 μm < T < 16.5 μm), and the theoretical threshold value for test 4 was defined as 16.1 μm. Since the precise threshold range was between 16.1 μm and 16.5 μm (16.1 μm < T < 16.5 μm), the scale of overestimation in practice was negligible. Hence, the revised upper limit of 16.5 μm was used as a proposed threshold for regular low-density samples. For the corresponding threshold model (TM3) modified from TM1 and TM2, the x in Eq. (7) was taken to (D + 3h)/4.
The outcomes of COMSOL simulation were validated by actual experiments. Specifically, we evaluated TM3 with 16.5 μm PS microspheres (Fig. 3d). The highest PRR was found at outlet 3 (77.6 ± 6.2%), while the rest of the microspheres was found at outlet 4 (22.4% ± 6.2%). Although the hPRR of 16.5 μm microspheres at outlet 3 were higher than that of 15 μm microspheres (52.4 ± 4.6%), it was lower than that of 21 μm microspheres (91.6 ± 7.3%) (Fig. 3e). Thus, TM1 remains the most accurate model for regular low-density particles in this study.
Evaluation of all models with irregular low-density particles. One of the main aims of redefining the model used to estimate the separation threshold is to expand the application range of inertial focusing techniques for samples of different densities, shapes, and sizes. Here, we investigated the accuracy of all three threshold models for irregularly shaped samples using polystyrene (iPS) particles of various shapes and sizes (< 100 μm). Through actual experiments, we obtained the recovery rates of 15 μm, 16.5 μm 21 μm iPS samples (± 2 μm) in outlets corresponding to the thresholds of TM1, TM2, and TM3. The hPRRs of 21 μm, 16.5 μm and 15 μm iPS particles were 41 ± 4.8% (outlet 3), 31.2 ± 10% (outlet 4) and 27.6 ± 4.9% (outlet 3) (Fig. 4). Then, we evaluated whether TM1 will underestimate the threshold by obtaining the PRR of the larger iPS particles Evaluation of all models with irregular high-density particles. The focus position of particles within the channels depends on the particle density 12 . However, since inertial microfluidics research is mainly designed for biomedical applications, low-density particles such as polystyrene microbeads ( ∼ 1.05 g/ml) were often used because their densities are similar to that of biological samples. In recent years, scientists are now interested in exploring the use of inertial microfluidics for higher density targets. Hence we used irregular talcum (iTalc) particles ( ∼ 2.75 g/ml) to investigate how particle density could affect the accuracy of the existing threshold models. It could be seen from the results that the hPRRs of 21 μm, 15 μm and 16.5 μm iTalc particles at outlet 3 were 38 ± 11.5%, 31.3 ± 9.1%, and 32.2 ± 9.9%, respectively (Fig. 5). These low PRRs of respective diameters were probably due to the high density of talcum. According to Eqs. (3) and (6), F L and F D are independent of the change in particle density, and thus an increase of particle density might not influence the balance of F L and F D . However, in inertial microfluidics, apart from the F L and F D , the particles migrating in the curved channel were also subjected to the centrifugal force F c which is given by: where ρ p and ρ f are the particle and fluid densities, respectively, V is the tangential velocity of the particles, and r is the radius of the particle trajectory 10 . In most studies, F c is usually ignored for proximate densities of particle and fluid medium ( ρ p − ρ f ) 12 . For regular and non-regular PS samples, since the densities of PS ( ∼ 1.05 g/ml) and water ( ∼ 1 g/ml) are very similar, F c acting on microbeads is negligible, and thus particle focusing depends on the balance of F L and F D (F L = F D ). However, for higher density particles such as talcum, the F c acting on the particles would be higher (keeping other parameters unchanged), and thus F c became non-negligible. When the particles were subjected to Dean vortices in the inner cross-section, F c acted in the direction opposite to F L , and therefore a larger F L was required to balance the additional F c and to focus the particles of higher density (F L = F D + F c ). Hence, the deviation of PRRs values from TM1 estimation strongly suggested the need for a more precise model for high-density targets.
Modification of TM1 for irregular high-density particles. Since larger particles will experience a larger F L in the system, the threshold for iTalc particles was expected to be larger compared to that of iPS particles. In order to derive a more precise model for samples with varying densities, we obtained the PRR of the larger iTalc particles (26 ± 2 μm) at each outlet. The results showed that the highest PRR of 26 ± 2 μm iTalc www.nature.com/scientificreports/ particles (50.2 ± 12.6%) found at outlet 3 was significantly higher than that of 21 μm (38 ± 11.5%) (Fig. 5d). This result was in line with our expectation that the larger F L induced by larger particles would balance the additional F c . Then, by analyzing the PRR of iTalc particles larger than 26 ± 2 μm (28-50 μm), we assessed whether the threshold was saturated at 26 ± 2 μm. The hPRR for larger iTalc particles at outlet 3 was found to be 48.6 ± 15.3%, which was very similar to the PRR of the smaller iTalc particles (26 ± 2 μm) (Fig. 5e). Therefore, the threshold for high-density particles (iTalc) was probably saturated around 26 μm. Similar to the result of iPS particles, a slight increase in size did not increase the highest PRR, despite the larger F L experienced by larger particles. This could be due to the presence of other opposite forces, F c and F D, which would also be increased in the system (F c ∝ a 3 and F D ∝ a ). Therefore, there was a saturation point at which the focusing effect did not increase as the particle dimension increased. To allow estimation of PRRs for samples with varying densities, we thus derived a model (TM4) as: where c was a density coefficient greater than or equal to unity (c ≥ 1). The value of c increased with the density of sample, for example, for polystyrene particles ( c = 1; ∼ 1.05 g/ml) and for talcum particles ( c = 1.215; ∼ 2.75 g/ml). Since the c value is equal to 1 for PS, and 1.215 for Talc, substituting above c values to TM4, the thresholds for iPS and iTalc particles were calculated to be 21 μm and 26 μm, corresponding to the highest PRRs obtained previously (Figs. 4 and 5).
Applicability of TM4 for irregular particles of different densities. In order to study whether TM4 could be applied to particles of different densities, we hypothesized that the density coefficient c of TM4 would increase linearly with particle density, from 1.05 g/ml to 2.75 g/ml (Fig. 6a), whereby the proposed linear correlation can be modelled by: where c is the density coefficient, and ρ is the particle density, respectively. We aimed to verify this linear correlation using irregular polyethylene terephthalate (iPET) ( ∼ 1.38 g/ml) and polyamide (iPA) ( ∼ 1.15 g/ml). According to Eq. (11), the density coefficient (c) of iPET was around 1.04. Based on TM4, the threshold for iPET was approximately 22.5 μm. We then obtained the PRRs of iPET particles of various sizes (22.5 ± 2 μm, 26.5 ± 2 μm, 20.5 ± 2 μm, and 28.5-50 μm). We found that the highest recovery rates of 22.5 μm, 20.5 μm, and 26.5 μm iPET particles were around 46.9 ± 2.3%, 42.1 ± 12%, 43.9 ± 4.3% and 46.2 ± 5.1% respectively (Fig. 6b-e). As the highest PRR of 22.5 μm particles was higher than that of 20.5 μm particles and close to 26.5 μm and 28.5-50 μm iPA particles, 22.5 μm was determined as the threshold for iPET particles. Then, iPA particles were  (Fig. 7a-d) showed that the threshold for iPA was around 21.5 μm (60% hPRR at outlet 3) because the hPRR of 21.5 μm iPA was higher than that of 17.5 ± 2 μm and similar to that of 25.5 ± 2 μm. Therefore, the linearity was validated (Fig. 7e).
Overall, we derived an improved model (TM4) to estimate the separation efficiency of samples with varying densities by incorporating the factor of particle density. Our studies demonstrated that TM1 was only suitable for the low-density particle (1.05 g/ml), while TM4 was applicable for samples with a broader range of densities, including high-density samples (2.75 g/ml). The linear correlation between c and ρ makes TM4 more applicable to samples of different densities. For example, since the value of c is 1 in TM4 for PS particles, TM4 will be exactly equivalent to TM1, which was the most accurate for PS particles, and thus TM4 retained the same accuracy as compared to TM1 for low-density samples. For intermediate-density and high-density particles (e.g., PET and Talc), since the centrifugal force acting on them was non-negligible, the parameter was considered in the derivation of TM4 by introducing a density coefficient c, which was customized for particles of different densities. Specifically, the threshold of suitable density range for TM4 was between 1.05 to 2.75 g/ml. We hope these studies can also provide a systematic approach for future efforts to derive an optimal prediction model for a wider range of biological samples.

Discussion
Some parameters essential to channel design, such as C L and U D in Eq. (3) and (5) can only be estimated through experimental data 12 . These equations can be modified by utilizing the size-dependent nature of R f (R f = F L /F D and R f ∝ a 3 ). However, there is a lack of standardisation of the threshold model used, leading to a significant discrepancy in threshold prediction. In addition, current models were developed using low-density samples ( ∼ 1.05 g/ml) and may not be applicable to samples of higher density. We believe that the applications of inertial microfluidics will be extended to the use of high-density particles in the foreseeable future. One potential application is to facilitate the monitoring of heavy metals in drinking water. These high-density metals (such as Ag, Cu, Ni, and Cr) are carcinogens present in very low concentrations. Nonetheless, current point-of-care technologies for monitoring metals are subjected to high limits of detection (LOD), of which the LOD values are higher than the maximum allowable limit suggested by WHO 24 . Therefore, an improved threshold model Figure 6. The proposed linear correlation between density coefficient, particle density, and the recovery rate of irregular polyethylene terephthalate (iPET) of various diameters. (a) The graph illustrates that c linearly increased with ρ within the range of 1.05 g/ml to 2.75 g/ml. The PRRs of (b) 22.5 ± 2 μm, (c) 20.5 ± 2 μm, (d) 26.5 ± 2 μm, and (e) 28.5-50 μm were indicated respectively. Since the hPRR of 22.5 μm was similar to 26.5 μm but higher than that of 20.5 μm, 22.5 μm was determined as the threshold for iPET. Data were shown as mean ± SD of three independent experiments. *p < 0.01, **p < 0.001, ***p < 0.0001. Scientific Reports | (2020) 10:18984 | https://doi.org/10.1038/s41598-020-75878-w www.nature.com/scientificreports/ incorporating sample density for point-of-care inertial microfluidic biosensors can be useful for the detection or removal of toxic metals from drinking water. Overestimation of the threshold, even as small as a few microns, will also significantly affect the design of a microfluidic device and prevent the optimal separation of the subpopulations of cells. For example, the diameter of white blood cells ( ∼ 9.5 μm) is similar to that of red blood cells ( ∼ 7 μm). Here, we evaluated the accuracy of the existing models (TM1 and TM2) and derived an improved model (TM4) to estimate the separation efficiency of samples with varying densities by incorporating the factor of particle density. TM4 provides a more precise estimation of the separation efficiency of samples with varied densities and will aid researchers in their design of microchannels for particle separation based on the inertial focusing effect.
Similar to TM1 and TM2, it should be highlighted that TM4 is only applicable to channels with a rectangular cross-section. Interestingly, the previous threshold models did not consider geometric considerations, and both TM1 and TM2 were derived based on rectangular cross-section channels. Since the focus of this manuscript is to expand on existing threshold models (TM1 and TM2) and incorporate the component of sample density for inertial microfluidics, the standard rectangular cross-section design was used for our simulation and experimentations. Our preliminary studies (Fig. S2) suggested that all current threshold models (including TM4) were indeed not applicable for microfluidic channels with varied geometric designs (e.g., trapezoidal). Future efforts should be directed towards the modification of existing models for different geometric shapes.
Other factors can also influence the sorting efficiency of particles. Previous studies have also emphasized that unless the asymmetry is large, the shape of the particle has only a small effect on inertial focusing, although the relevant mechanism was not clarified in the literature 25 . In order to show whether there were some PS powders with great asymmetry, a scatter diagram of the circularity distribution of iPS powders was generated (Fig. S3). Since the scattering points were uniformly distributed from 0.1 to 1, the shape of the particles was highly diverse. This implied that a certain portion of particles could have a higher degree of asymmetry. As previous research on the correlation between the particle shape and the focus position was carried out with straight channels, further studies will be required to consider how Dean resistance is affected by the shape of the particles.
The experimental results of the microbeads were different from that of the simulated runs. For example, experimental results showed that outlet 4 had a higher PRR of 21 μm PS microbeads, which could be due to the stronger F L experienced. On the contrary, the 16.5 μm microbeads were less focused in outlet 3, possibly due to a moderate F L acting on the microbeads 26 . The shift of the equilibrium position from outlet 4 to 3 in experimental results could also be due to differential positioning depending on the nature of F L . Particles circulating along with Dean vortices at the inner section would experience larger F L at the proximate center location than at proximate inner channel wall 7 . Compared with 21 μm microbeads, 15 μm and 16.5 μm microbeads could experience Figure 7. Recovery rate of irregular polyamide (iPA) particles of (a) 17.5 ± 2 μm, (b) 21.5 ± 2 μm, (c) 25.5 ± 2 μm, (d) 27.5-50 μm diameters. As the highest recovery rate of 21.5 ± 2 μm iPA (59.8%) was higher than that of 17.5 ± 2 μm (52.1%), and very close to that of 25.5 ± 2 μm (55.7%), 21.5 μm was determined as the threshold for iPA. Data were shown as mean ± SD of three independent experiments. *p < 0.01, **p < 0.001, ***p < 0.0001. (e) The linearity of density coefficient and particle density. Points A, B, C, and D corresponded to the experimental data obtained with iPS, iPA, iPET, and iTalc particles, respectively.

Scientific Reports
| (2020) 10:18984 | https://doi.org/10.1038/s41598-020-75878-w www.nature.com/scientificreports/ sufficient F L to balance F D when they were located in a lateral position closer to the center of the channel. In this study, COMSOL was only a means to provide a rough estimation, as the accuracy of COMSOL was also influenced by many factors, such as the number of mesh elements. Increasing the number of elements generally led to higher accuracy, but the demand for computation power was also higher. In this study, the number of mesh elements was standardised to around 37,000. Despite these limitations, we demonstrated that the outcomes of the numerical simulation were close to our actual experiments of 15 μm and 21 μm microspheres. Biological samples are often heterogeneous, especially when clinical samples are involved. Prior studies on the detection of mesenchymal cancer cells from bladder wash urine and blood samples had demonstrated the fluctuations in target cell retrieval in clinical samples 23 . Here, we evaluated the performance of different models using particles of different size, shape, and density. Our result demonstrated that TM1 was a more accurate model for low-density particles. However, PPRs of higher density samples deviated from expected values. We thus modified TM1 by considering the effect of the additional centrifugal force F c acting on high-density particles in the curvilinear channel. The derived model (TM4) comprises of a density coefficient c such that the threshold obtained fits experimental values. The robustness of TM4 for samples of varying densities (1.05 g/ml ≤ ρ ≤ 2.75 g/ml) was confirmed with intermediate-density particles (PET and PA). In addition, TM4 was also as accurate as TM1 for low-density particles (Table S3).
Future analysis should cover other parameters that could influence the threshold, such as changes in channel aspect ratio, channel radius (R), and cross-section dimensions (same aspect ratio), which were reported to affect inertial focusing 27 . For example, apart from shape and density, particle deformability is one of the factors affecting the inertial effect owing to an additional life force experienced by non-rigid particles 28 . Similar to other inertial studies using polymeric particles (usually PS), all particles were assumed to be rigid because of their relatively high elasticity ( ∼ 2 − 4GPa ). However, the mechanical characterization of cells is similar to a drop of liquid enclosed by a membrane exhibiting low elasticity 29 . For example, the elastic modulus of mammary cells (MCF10A) and breast cancer cells (MCF7) widely used for studying breast cancer were around 440 Pa and 300 Pa, respectively 30 , while the stiffness of healthy leukocytes was even lower (10-60 Pa) 31 . Based on differences in deformability, isolation of more deformable MCF7 and healthy leukocytes from less deformable MCF10A and immature leukocytes had been reported in previous studies despite their similarity in sizes 31,32 .
Overall, we demonstrated that TM1 remained accurate for regular and non-regular low-density particles (e.g., PS). However, a modified model TM4, which incorporated the density coefficient c (TM4), was more suitable for samples with varied densities, especially those with high-density particles (e.g., Talc). The linear correlation of density coefficient c and particle density ρ was also validated. These findings are highly applicable to applications where the density of targeted particles is high or heterogeneous.

Methods
Device feature and fabrication. The spiral device with a rectangular cross-section (the main device in this study) consisted of a 10-loop curvilinear channel with 300 μm spacing between loop interval. The width and height of the channel were 500 μm and 220 μm, respectively. The device had 1 inlet and 5 outlets. The end of the 500 μm wide channel was evenly split into five outlets, each of which had 100 μm width. Spiral devices with trapezoidal cross-section consisted of an 8-loop curvilinear channel with 300 μm spacing between loop interval. The width of the channel was 600 μm, while the inner and outer heights were 80 μm and 130 μm, respectively. The device had 1 inlet and 2 outlets. The 600 μm wide channel was evenly split into two outlets.
Microdevices were fabricated using standard soft lithography from an aluminum master mold. Briefly, polydimethylsiloxane (PDMS) was prepared by mixing a 10:1 ratio of base and curing agent (Dow, Germany). The prepared PDMS was degassed in degasser and poured onto the master mold. The mold filled with prepared PDMS was baked at 70 °C oven for 2 h. After baking, the solidified PDMS was peeled off from the mold carefully. The locations for inlet and outlets were punched using a 1.5 mm diameter puncher (Integra, USA). The punched PDMS device was bonded to a flat solidified PDMS layer to complete the channels in a plasma machine. After plasma treatment for 3 min, the bonded PDMS device was placed in a 70 °C oven for 30 min to strengthen the bonding. Finally, 15 cm and 6 cm long sections of flexible plastic tubing (Tygon, USA) were inserted into the device inlet and outlets, respectively.
COMSOL flow simulation. To evaluate the particle focusing of different sizes for respective threshold models. COMSOL Multiphysics software was used to simulate the inertial focusing on our existing microdevice. The three-dimensional (3D) image of the microdevice was created by AutoCAD in STL format (Fig. 1a). All dimensions were the same, except for the shorter length of the branch at all outlets (to save computation time). The properties of water were applied to systematic fluid for simulation (ρ = 1000 kgm −3 ; µ = 10 −3 kgm −1 s −1 ). The laminar flow (spf) module with incompressible flow property was used in the system. F low velocity was set at the inlet, while the fixed pressure boundary condition was set at all outlets. Using these perimeters, a flow magnitude profile for the system was simulated (Fig. 1b).
Based on the simulation of the flow profile, particle tracing for the fluid flow (fpt) module was chosen to analyse the motions of particle trajectories. The drag force, wall-induced life force, and shear gradient life force parameters were applied to the system. Interactions between particles were ignored as the diluted sample was often used in the inertial microfluidic application. Solid particles ( 1.05 g/ml) was used to model spherical polystyrene (PS) microbeads. Each set of samples consisted of 3 particles of the same size evenly distributed at the inlet to initiate simulation (Fig. 1a). The simulation time for each set of particles was 2 s with 0.01 s time step for tracing particle trajectories. The 'freeze' condition was applied at the end of outlets to ensure all particles within the boundaries. The number of particles freezing at the end of boundaries was counted.
Scientific Reports | (2020) 10:18984 | https://doi.org/10.1038/s41598-020-75878-w www.nature.com/scientificreports/ Flow rate optimization. For the inertial microfluidic system, the Reynolds number (Re) is required to be within the range of 1 to 100 (~ 1 < Re < ~ 100), which corresponds to the volumetric flow rates approximately from 0.021 ml/min to 2.1 ml/min. Since a higher flow rate was preferable, different flow rates (1.5 ml/min, 1.7 ml/min, 1.9 ml/min, and 2.1 ml/min) were used to process 21 μm PS regular microbeads in the system. 21 μm spherical microbeads were used as the particle size was above the thresholds of both TM1 and TM2, which indicated that 21 μm PS microbeads should be well-focused in the device. We confirmed that based on the highest PRR (91.7% PRR) obtained under this parameter, the optimal flow rate is 1.7 ml/min (Fig. S4). Although the PRR derived from the flow rate of 1.9 ml/min suggested that this parameter could also be used, it was close to the limit of Reynold's number applicable to inertial focusing based microfluidics, and thus only 1.7 ml/min was used for subsequent experiments.
Preparation of irregular particles. iTalc (baby powder) was commercially available (Johnson & Johnson China Ltd, China). For iPS, iPET, and iPA powders, they were prepared using a Retsch CryoMill cryogenic grinder (Haan, Germany). After a 7-min precooling stage with five shakes per second, the grinding stage was used for 1.5 min with 25 shakes per second in the presence of liquid nitrogen. The crushed microplastic powders (< 100 μm) were then collected through a 100 μm filter. The powders were mixed with DI water for experimentation in a biological safety cabinet to minimize potential contamination.
Particle sorting experiments. 1 μl microbead suspension was mixed well with 1 ml deionized water and transferred to 10 cc syringe. Another 10 cc syringe was prepared by injecting 2 ml deionized water and was loaded to a syringe pump. The inlet and outlets of the microdevice were connected to the syringe and five centrifugal tubes using plastic tubing. The device and centrifugal tubes were arranged accordingly (Fig. 1c). The syringe pump was run at 1.7 ml/min for 10 s. The loaded syringe was replaced by another syringe with microbead suspension. The syringe pump (New Era Pump Systems Inc., USA) was run at 1.7 ml/min for 20 s. The sample collected at centrifugal tubes was then transferred to a 24-well plate (SPL Life Science, Korea). Each well with microbead suspension was mixed well with an additional 100 μl deionized water. The 24-well plate was placed at rest for 30 min to settle down the particle suspension. The number of regular particles in each well was counted under an optical microscope (Nikon, Japan) (Fig. S5). All experiments were carried out in triplicates.
For irregular particles, all images captured by the optical microscope (Fig. S6) were analysed by ImageJ. Briefly, images were converted to 8 bits with suitable threshold levels to customize the particle counting function. The adjusted images were carefully compared to their original images to ensure that the threshold selected was optimal. Their circular diameters were measured with ImageJ and Excel. Briefly, we first obtained the boundary area (A) of each particle by ImageJ. Then, the circular diameter (d) was calculated by d = √ (4*A/π) using Excel. Due to the high variation of sample sizes, the upper and lower ranges of two microns (± 2 μm) were included for the measurement of the specific dimension.


extfieldsofstudy:
Physics
Computer Science

text:
A survey on key technologies of privacy leakage detection for Android platform

With the rapid advancement of the mobile internet, smartphones have been the essential tools for communication and entertainment. In recent years, Android, developing quickly, has been the largest mobile operating system in the world. However, due to mass private data generated by Android intelligent devices, privacy disclosure is a serious problem in security for the Android platform. Based on the relevant background, this thesis indicates the problems of privacy protection that the Android platform faces, summarizes the detecting techniques for privacy disclosure on the Android platform at the present stage, and introduces the research status and progress with the aim of offering a reference for relevant practitioners and researchers.

The Outstanding problem of privacy protection on the Android platform
On the aspect of privacy protection of users, the Android platform takes a series of security measures to protect the privacy security of users. In order to protect smartphones from losing or being stolen, android provides the lock screen and full disk encryption for data partition of the Android system by using master keys which are randomly generated. Therefore, if users lose their smartphones, the privacy information storing in the smartphones cannot be disclosed unless the person involved knows the password of the lock screen. As for the Apps operating on the Android platform, Android system also provides a series of security measures such as the sandbox mechanism, the permission mechanism, and the signature mechanism. In terms of the permission mechanism, Android system, based on the permission management, provides a set of privacy protection mechanism to protect private data for users. In the privacy protection mechanism, App developers must make a clear statement of the required permissions in the document AndroidManifest.xml and strictly follow the security specification set by Android system to avoid these permissions being disclosed to the malware. Before installing the new App, users should carefully read all permission statements of the App, and only by being accepted by users, the App can be installed in the devices of users. Android system should make sure that only with the relevant permissions, Apps can be allowed to have access to private data or related functions. However, common users do not know whether the permissions applied by Apps are safe or not. As for Apps, Android system also provides the protection mechanism of digital signature. When users install an App, the system can detect the digital signature of the App. The system will detect whether the App was once installed in the smartphone or not, which means whether the name of the program package about to install is the same as that of the installed program package or not. If the names are matched, the system will treat the App about to install as the upgrade version of the installed App and then replace it. If the names are unmatched, the system will stop installing. Although some protection mechanisms such as sandbox mechanism and SLINUX mechanism can protect users' privacy to some extent and restrain the rampancy of the malware, these measures show the weakness when they tackle endless black technologies and attack techniques. Under this circumstances, in the study of how to improve the privacy protection of Android system, there are some outstanding problems to solve.
(1) Privacy data hiding and its deniability. Because the Android system only protects the inner privacy data by password of lock screen, and most people adopt easy and simple patterns to lock the screen, when users' smartphones are lost and users are hijacked, it is easy for others to get the passwords which may result in the inner privacy disclosure. If the system can encrypt to hide the private data in the smartphones, it can fill the gap of the weak protection of passwords of lock screen. By this way, even the password is cracked, or users are forced to give the password, the hidden private data can be well protected and users can deny the private data.
(2) The safe and reliable system operation environment. Currently, most of Apps on the Android platform contain service components, running in the background. This will cause a large number of unreliable Apps (such as entertainment or game software) operating when users use sensitive Apps like Alipay. Much worse, there are some malware gaining root privilege by making use of system bugs. Under this condition, sensitive Apps will be attacked by the malware so that the account information will be disclosed. Additionally, the private data storing in the smartphone may be stolen by the malware.
(3) Permission control of fine grit of internal storage. At present, with the price decline of the flash memory, the capacity of internal memory of smartphones rapidly increases to meet the needs of user information storage, which makes device manufacturers gradually no longer to support the removable TF memory. As for users, their personal data can only be stored in the internal memory. App developers usually store some key data in the memory (emulation memory). However, in order to have access to the internal memory and all files in the emulation memory, Apps just need to apply for the permissions of SDCARD_R and SDCARD_W. This is the main way of privacy disclosure, also the main target attacked by the ransomware with mobile encryption.

The Static Detecting Technique for Privacy Disclosure
As the common analysis technique in the field of mobile security, the static analysis technique for detection in Android means that instead of a code of an App under detection being operated, some techniques like lexical analysis and syntax analysis are used to analyze and detect the source code, decompiled code or intermediate code of the App to make sure whether the App will give rise to the privacy disclosure or not. The static detecting technique stresses on static state, so the main way during the process of detection is analyzing and scanning the decompiled code. There are two situations in the static detecting technique of Android Apps. The first one is to analyze and detect the Dalvik bytecode files or Smali files made by decompilation. Another is to analyze and detect the Java code made by decompilation. The key to distinguish these two situations is the different objects of analysis. The static analysis technique can be divided into two analysis methods, one of which is based on the feature code while the other is based on the syntax and semantic of the code.
(1) The analysis method based on the feature code is to generate the feature code for malicious act by analyzing the malicious App which is known to all. In fact, analyzing the Apps which were not analyzed before means to match the feature codes. Compared with the dynamic analysis technique, it is unnecessary for the static analysis technique to use the extra resources such as sandbox and virtual machine, and before the operation of malware, it can detect and scan. For this reason, it is characterized by high efficiency and low risk. However, the result of this kind of analysis method relies on the completeness of the feature data base. With the new variation viruses appearing constantly, the feature data base will be larger, which will not only lead to the detecting efficiency decline, but also increase the expense on maintaining the enlarged feature data base.
(2) The analysis method of code semantic is the detecting method based on the syntax and semantic of the code, which scans the code of malware and acquires its implementation as well as the semantic relation according to the features of Java syntax and semantic to avoid some interference like code obfuscation. The basic process of this kind of detecting method is to analyze the framework of APP codes and walware codes, and then extract the syntax and semantic features of both codes. After that, the semantic association diagrams will be set up based on these features and the similarity of two diagrams will be calculated.
To sum up, the general process of the static analysis for Android Apps is shown as Figure 1. It can be seen that the static analysis method is to scan and analyze the control flow and data flow to ensure whether it will lead to the privacy disclosure or not when the Android Apps are not under operation. AndroidLeaks However, these analyzing and testing systems mainly aim at the scanning and analysis of Java source code. At present, most of Android Apps are not open source. Therefore, firstly, there are many measures like decompiled operation for Apps and reverse engineering should be taken to extract the features of the Apps and get the information like relation of function call and binary sequence. In literature [6], FlowDroid provides more accurate information for the transformation of life cycle of Apps by carefully analyzing the whole life cycle of the Apps. Through developing DroidBench on the Android platform, FlowDroid can test the accuracy of taint detecting methods. The result shows that its accuracy is higher than AppScan and Fortify. The sandbox mechanism serves as a mechanism running in isolation. When it detects the privacy disclosure of Apps, it needs to detect the privacy leakage source of the system and API of the privacy leakage point in real time. Once the App invokes these interfaces, records should be taken in real time. Because the sandbox mechanism needs additional resources and the monitoring points are discontinuous, there lacks logical relationship and the context. In the analysis of privacy disclosure made by the dynamic taint tracer technique, users' sensitive information are often marked as the tainted data and then the spread track will be traced. If the tainted data are found to be leaked by message, the internet or bluetooth, the App is judged to exist the privacy leakage problem. The detailed process is shown as Figure 2.  [7], the privacy disclosure testing tool developed by Enck et al. uses the dynamic taint tracer technique early on the Android platform. According to TraintDroid, it firstly marks some privacy leakage sources (such as short messages, geographical location information and IMEI), mainly adding the 32-bit tainted label near the location of each data variable. Secondly, it traces the taints of coarsness (the level of IPC and file messages) and fine grit (the level of Dalvik virtual machine instructions). The last step is to detect API at the exit of the internet, files and short messages, including detecting the private messages about to send to judge whether they are the tainted data, recording tainted private data and the label type, and removing the tainted label before sending messages. However, the disadvantage of TaintDroid is that the process of the analysis detection requires human intervention rather than automatically detect privacy disclosure.
At present, most of safe-defend software adopt the dynamic detecting technique [7]. For example, in order to detect privacy disclosure, Paranoid Android [8], developed by Georios et.al, copies the operation condition of Android Apps, checks and kills the virus on the remote server, and detects the taint spread by taking advantages of the functions of the virtual machine including recording and playback. Homyack et.al upgrades TaintDroid as the tool which can filter the private information when the information is disclosed. It can largely protect the users' private data without affecting the operation. The literature [10] [11] belowdesigns and finishes the sandbox system on the Android platform based on the detection of basic dynamic act. Additionally, it can also monitor the operating act of Apps in the frameworks of different levels. The DroidScope system in literature [12] makes a dynamic detection for malware on the Android platform by recreating SalvikVM of the Android system.

Conclusions
Combined with the relevant background of the Android platform, the thesis points out the outstanding problems about privacy protection which needs to be solved by the Android platform, and summarized the research status based on the private leakage detecting techniques on the Android platform, including the static detecting technique and the dynamic detecting technique. In General, various privacy leakage detecting techniques have been the focus of researchers, and they needs to improve to fill the gaps even though they grow rapidly. With the fast execution and the comprehensive coverage of codes, the method of static analysis yet cannot solve the problems like code obfuscation and high rate of false alarm while 2019 2nd International Symposium on Big Data and Applied Statistics Journal of Physics: Conference Series 1437 (2020) 012006 IOP Publishing doi:10.1088/1742-6596/1437/1/012006 5 the method of dynamic analysis can better solve the problems of the static analysis, but it also exists some problems like insufficient coverage rate of codes and low operation. In conclusion, how to combine the advantages of the static analysis with the dynamic analysis to make a mixed analysis is the trend of privacy protection on the Android platform in the future.