Abstract

Releasing evolving networks which contain sensitive information could compromise individual privacy. In this paper, we study the problem of releasing evolving networks under differential privacy. We explore the possibility of designing a differentially private evolving networks releasing algorithm. We found that the majority of traditional methods provide a snapshot of the networks under differential privacy over a brief period of time. As the network structure only changes in local part, the amount of required noise entirely is large and it leads to an inefficient utility. To this end, we propose GHRG-DP, a novel differentially private evolving networks releasing algorithm which reduces the noise scale and achieves high data utility. In the GHRG-DP algorithm, we learn the online connection probabilities between vertices in the evolving networks by generalized hierarchical random graph (GHRG) model. To fit the dynamic environment, a dendrogram structure adjusting method in local areas is proposed to reduce the noise scale in the whole period of time. Moreover, to avoid the unhelpful outcome of the connection probabilities, a Bayesian noisy probabilities calculating method is proposed. Through formal privacy analysis, we show that the GHRG-DP algorithm is -differentially private. Experiments on real evolving network datasets illustrate that GHRG-DP algorithm can privately release evolving networks with high accuracy.

1. Introduction

As more and more individual information in social network is published, releasing network data might pose threats to individual’s privacy. To protect the privacy of individual information, network data should be sanitized before releasing. Differential privacy has been proposed to address such problem. Unlike the anonymization-based private algorithms (e.g., k-anonymity [1] and l-diversity [2]), differentially private algorithms are randomized algorithms. Differentially private networks releasing algorithms protect the essential structural information and ensure that the connection relationship of each participant is insensitive.

In the differentially private static network literature, the standard technique is to add noise to the output of an otherwise non-private algorithm. However, small changes of the network structure should cause a lot of impact point to query answers and it leads to a large amount of perturbation noise. Thus, synthetic network generation methods are used to substitute these standard methods. Roughly speaking, such synthetic network generation methods involve two steps: (1) using the real network to learn the parameters for a model of the network; (2) using the learned model to draw synthetic network from its probability distribution. These synthetic network generation methods could “dilute” the impact of small changes of the network structure by capturing the connection probabilities between vertices.

However, the interactions between vertices are often dynamic in nature, and how the network changes deserves to be taken into account. Given a sequence of evolving networks, the standard technique of differential privacy is to add noise to each snapshot of the evolving networks. In this way, each snapshot of the evolving networks satisfies differential privacy individually and turns the static differential privacy algorithm into a “dynamic” version. However, this kind of method suffers from poor performance. The root cause is it leads to a large amount of repetitive perturbation noise. We observe, in practice, the network structure between sequent snapshot changes within a local area. The remaining part of the network is almost unchanged and adding repetitive noise to the same connection probabilities is inefficient. Furthermore, the connection probabilities between vertices of the network may obtain unhelpful outcomes when the connection probabilities equal 0 or 1.

To this end, we propose a novel differentially private evolving networks releasing algorithm, called GHRG-DP (i.e., differentially private based on generalized hierarchical random graph). In particular, we first use generalized hierarchical random graph (GHRG) model to represent the structure of evolving networks and utilize a hierarchical structure (called dendrogram) of the GHRG model to record connection probabilities between any pair of vertices which have appeared over a period of time. Unlike the popular hierarchical random graph (HRG) model, GHRG allows tree nodes of dendrogram to have any number of children and quantify the change of network’s structure. As the network structure between sequent snapshots changes within a local area, it is important to design a method to generate a good dendrogram of each snapshot of evolving networks. We do not directly generate the best-fitting GHRG by MCMC method, but we generate the GHRG by adjusting the GHRG of the prior time. Then, we model each connection probability as a distribution instead of a point value and add noise to the parameters for the model of the network instead of adding noise to the point value of the connection probability. Finally, we compute connection probabilities’ posterior by Bayesian theory.

To make the GHRG-DP satisfy differential privacy, we first sample a GHRG by a Markov chain Monte Carlo (MCMC) method by exponential mechanism while satisfying differential privacy. Then, we adjust the GHRG by exponential mechanism while satisfying differential privacy as well. Thirdly, we add Laplace noise to the parameters for the model of the network. Through formal privacy analysis, we prove that GHRG-DP satisfies -differential privacy. In summary, we present several contributions:(1)We introduce generalized hierarchical random graph model as an effective means of encoding the evolving networks and propose a method to adjust the GHRG model which is suitable for evolving networks.(2)We develop a novel differentially private evolving networks releasing algorithm, GHRG-DP. To ensure that the releasing evolving networks under differential privacy do not incur excessive noise, we propose a method by adding noise to the parameters for the model of the network.(3)Through privacy analysis and experiments on real evolving network datasets, we prove that GHRG-DP algorithm can privately release evolving networks with high accuracy.

The rest of our paper is organized as follows. Section 2 provides a literature review. Section 3 introduces necessary background on differential privacy and hierarchical random graph model. Section 4 identifies the weakness of applying differentially private static network releasing algorithm to evolving networks. In Section 5, we describe our GHRG-DP algorithm in detail. Section 6 gives the privacy analysis. Section 7 reports the experimental results. Section 8 concludes the paper.

Many existing works about social network differential privacy focus on social network analysis. These methods output some network statistics under differential privacy such as degree distribution, subgraph number, and clustering coefficient. In particular, Dwork et al. [3] proposed a differentially private method to answer the queries by adding noise to outcome directly. Hay et al. [4] proposed a differentially private method in a postprocessing phase and computed the consistent input most likely to have produced the noisy output. Hay et al. [5] used this method to estimate the degree distribution. Nissim et al. [6] introduced the concept of smooth sensitivity to protect individuals’ privacy by adding a small amount of random noise to the released statistics, such as triangle count of a network. Zhang et al. [7] analyzed the released statistics through a ladder function and reduced the sensitivity effectively. Cheng et al. [8] presented a two-phase differentially private frequent subgraph mining algorithm called DFG. In DFG, frequent subgraphs are privately identified in the first phase, and the noisy support of each identified frequent subgraph is calculated in the second phase. Ding et al. [9] published the triangle counts satisfying the node-differential privacy with the triangle count distribution and the cumulative distribution. Sun et al. [10] formulated a decentralized differential privacy scheme named DDP, which requires that each participant considers not only their own privacy but also that of their neighbors involved in their ELV. Wang et al. [11] studied the problem of differential privacy for weighted network through the probability model.

Differentially private social network releasing algorithm also draws attention. Sala et al. [12] introduced a differentially private graph model called Pygmalion. Pygmalion extracts a graph’s detailed structure into dK-graph, introduces noise into the resulting dataset, and generates a synthetic graph. Lu and Miklau [13] proposed algorithms for privately estimating the parameters of exponential random graph models (ERGMs) of a network. Based on the idea of stochastic Kronecker graph model, Mir and Wright [14] used maximum likelihood estimation to privately estimate the parameters. Xiao et al. [15] proposed a differentially private network publishing method, HRG-MCMC, witch computes an estimator of a given graph in the hierarchical random graph (HRG) model in a differentially private manner and samples possible HRG structures in the model space via Markov chain Monte Carlo (MCMC) witch satisfies the exponential mechanism. Qin et al. [16] investigated techniques to ensure local differential privacy of individuals while collecting structural information and generating representative synthetic social graphs. Chen et al. [17] presented a method for publishing differentially private synthetic attributed graphs, which is able to preserve the community structure of the original graph without sacrificing the ability to capture global structural properties.

Many existing works also focus on ensuring dynamic network data privacy. Bhagat et al. [18] analyzed the privacy risks of publishing multiple instances of the same network independently and proposed methods to perform group-based anonymization. Tai et al. [19] introduced the concept of k2-degree anonymity, which limits the probability of a vertex being re-identified to 1/k. Krzysztof Juszczyszy measured the transitions during network evolution by link prediction. Medforth and Wang [20] proposed an anonymization method to solve “degree-tail” attack. Macwan and Patel [21] proposed an anonymization approach to preserve the user identity from all the published time-series datasets of a social network. Yue et al. [22] proposed an efficient and adaptive general graph anonymization framework for incremental data publication. They proposed an anonymization process restart issue (APRI) and designed an activated function to determine whether an anonymization process should be restarted at a certain time in order to solve the problem of high time complexity and large information loss in the incremental data publication. Zhiyuli et al. [23] attempted to model the hierarchical and dynamic features of social networks by designing a damping-based sampling algorithm corresponding to a local search-based incremental learning algorithm, which can easily be extended to large-scale scenarios.

Most existing differentially private graph models do not fit evolving network. In this work, we develop a differentially private network publishing method for evolving network.

3. Background

3.1. Hierarchical Random Graph

Suppose the input network dataset is a simple undirected graph , where V is the set of vertices and is the set of edges. Claustet et al. [24] proposed the hierarchical random graph (HRG) model. HRG uses a hierarchical structure and connection probabilities between any pair of vertices to denote a graph G. The hierarchical structure is represented by a dendrogram T. The dendrogram is a binary tree which has n leaf nodes corresponding to the n vertices. Each internal node r represents a connection probability . In particular, , where and denote the numbers of leaves of left subtrees and right subtrees , respectively. is the number of edges between the leaves of and .

Given a graph G, we use the likelihood of the HRG to measure how much it matches G. The likelihood can be calculated as follows:where is the Gibbs–Shannon entropy function. A higher likelihood corresponds to a better representation of the network’s structure.

3.2. Generalized Hierarchical Random Graph

Peel and Clauset [25] introduced the generalized hierarchical random graph (GHRG) model in 2014. It also defines a couple of data structures in formal. Unlike the popular HRG model, the GHRG allows the nodes in the dendrogram to have two or more child nodes. It models community structure at all scales and provides accurate fits to social networks. In the dynamic evolving systems, the GHRG could improve the interpretability for varying a network’s structure.

3.3. Differential Privacy

Given a graph G, differential privacy ensures that the outputs are approximately the same even if any edge is arbitrarily added or deleted in the graph. Thus, the presence or absence of any edge has a negligible effect on the outputs. We define two graphs and to be neighbors if they satisfy , , and . -differential privacy is defined as follows.

Definition 1. (-differential privacy [3]). A randomized algorithm A is -differential privacy if for any two neighboring graphs and , and for any output ,Differential privacy is based on the concept of global sensitivity of a function f. It is used to measure the maximum change in the outputs of f when any edge in the graph is changed. The global sensitivity of f is defined as .
Differential privacy can be achieved by Laplace mechanism and exponential mechanism. The Laplace mechanism is mainly used for functions whose outputs are real values. Differential privacy can be achieved by adding properly noise drawn randomly from Laplace distribution to the true answer.

Theorem 1. (Laplace mechanism [3]). For any function with sensitivity , the algorithmsatisfies -differential privacy, where are i.i.d. Laplace variables with scale parameter .
The exponential mechanism is mainly used for functions whose outputs are not real numbers. The main idea is to sample the output data O from the output space according to the utility function u. The global sensitivity of u is .

Theorem 2. (exponential mechanism [26]). Given a graph G and a utility function , the arithmetic whose output is with probability proportional to satisfies -differential privacy.

Theorem 3. (sequential composition [27]). If each arithmetic provides -differential privacy, a sequence of over the same database D provides -differential privacy.

4. A Straightforward Approach

In this section, we present a straightforward approach to public evolving networks under differential privacy. The main idea of this straightforward approach is to add noise to every snapshot of the evolving network and let the graph of every timestamp satisfy differential privacy. In particular, we first identify the hierarchical random graph (HRG) that best fits . presents the snapshot of evolving networks at time t. And then, we generate from the identified HRG. Recall that an HRG consists of a dendrogram T and associated probabilities set . According to the HRG-MCMC method in [15], we can sample a good dendrogram T by means of Markov chain Monte Carlo (MCMC) procedure which satisfies exponential mechanism. For the associated probabilities set, we add noise to them directly to generate . At last, we generate from and dendrogram T. To satisfy differential privacy, the acceptance ratio iswhere is the global sensitivity of the utility function. It is worth noting that the generating process of timestamp are independent of each other.

4.1. Privacy Analysis of Straightforward Approach

Firstly, at every timestamp, MCMC procedure under exponential mechanism satisfies differential privacy, and the sensitivity is . When n is even , and when n is odd. Then, generating satisfies differential privacy. Hence, based on the sequential composition, the method in each timestamp satisfies differential privacy. When the time span is t, the straightforward approach satisfies differential privacy.

4.2. The Disadvantage of Straightforward Approach

However, the straightforward approach produces poor results. Firstly, as we observe, the structure of evolving network changes a little most of the time. The structure and the dendrogram of a snapshot are similar to a contiguous snapshot. The straightforward approach samples dendrogram T via MCMC procedure in each timestamp, respectively. It results in a lot of repetitive computation and drastically reduces the utility of the results. Then, the associated probability may lead to calculation impairments when we calculate the likelihood function. Consider the case where zero connections , or all connections ; we have or 1; thus the likelihood drops to 0 and results in an unhelpful outcome.

5. GHRG-DP Algorithm in Description

We propose a differentially private evolving networks releasing method, called GHRG-DP. This method is based on generalized hierarchical random graph and can solve the deficiency discussed in Section 4. Our solution is composed of two schemes: (1) we propose a dendrogram construction method for evolving network, called dendrogram construction based on adjustment (DCBA) (see Section 5.1). DCBA samples a dendrogram by MCMC process at t time. To achieve differential privacy, we employ the exponential mechanism during the MCMC process. Then, we propose a dendrogram structure adjustment method. We structure a dendrogram of t + 1 time by adjusting the structure of locally. (2) Given the graph and the dendrogram , we propose a noisy probability calculation method based on Bayesian theory (see Section 5.2). In this method, will not be a single numerical value, but a distribution. Thus, we could quantify the uncertainty of and avoid equal to 0 or 1. Then, we calculate the posteriori distribution of by Bayesian theory. To achieve differential privacy, we add noise to the parameters or hyperparameters to achieve the noisy probability . At last, we place edges to with probability and get a privatized synthetic graph.

5.1. Dendrogram Construction Method: DCBA

Given evolving graphs , we measure how plausible this dendrogram is to represent the graph by the logarithm of the likelihood. This log-likelihood can be computed by Bayes theorem. In fact, we find that the structures of adjacent dendrogram are similar. For example, in dataset AS-733 [28], the log-likelihood of dendrograms in adjacent time changes less than 0.05. Thus, we could get through adjusting locally and reduce the amount of noise we need to add.

To this end, we propose a locally dendrogram structure adjustment algorithm, called DCBA (dendrogram construction based on adjustment). The main idea of DCBA is structuring the dendrogram by GHRG model to satisfy the dynamic. When the structure of changes a little compared to , we could get through adjusting locally. Otherwise, we should sample a new through the MCMC algorithm which satisfies the differential privacy. In DCBA, we adjust the structure of dendrogram locally in order to simplify process of structuring the dendrogram. In particular, when a new vertex s is added to , based on the GHRG model, we could first, respectively, insert s into as the brother node of each leaf node and get a candidate set of dendrogram. For example, Figure 1 gives an example of the candidate set of dendrogram. Then, we compute the log-likelihood of each candidate. Finally, we choose the dendrogram which has the maximum log-likelihood value to be .

The DCBA algorithm is shown in Algorithm 1. Given the evolving graphs , we first judge the change situation of the graph structure at time t + 1 through change-point detection method (line 1). We use the method from [25] to detect the change point and choose threshold. Next, when the change of the graph structure is greater than the threshold, we should sample the dendrogram renewedly through MCMC which simulates the exponential mechanism (lines 2–10). Let the utility function be , where is the sensitivity of the utility function. On the contrary, when the change of the graph structure is less than the threshold, we should adjust the structure of dendrogram locally. We get a candidate set of dendrogram and compute the log-likelihood of each candidate (lines 12–15). Finally, we choose the dendrogram through the exponential mechanism (line 17). Let the utility function be and we sample with probability proportional to . is the entire candidate set of the new dendrogram. It is easy to know that the global sensitivity and we set it as 1.

Input:sampled dendrogram at time t, evolving graphs , privacy parameter , .
Output:sampled dendrogram at time t + 1.
(1)if the result of change-point detection of is true then
(2) Initialize the Markov chain by choosing a random starting dendrogram ;
(3)  for Markov chain step do
(4)  Randomly pick an internal node r in ;
(5)  Pick a neighboring dendrogram of by randomly drawing a configuration of r’s subtrees;
(6)  Accept the transition and set with probability ;
(7)end for
(8)if equilibrium is reached then
(9)  Return the sampled dendrogram ;
(10)end if
(11)else
(12)for each new node do
(13)  for each internal node u whose child nodes are leaves in do
(14)   regard s as u’child node and construct a candidate set of the new dendrogram;
(15)   compute the log-likelihood of each candidate ;
(16)  end for
(17)  sample the new dendrogram with probability proportional to from the candidate set;
(18)end for
(19)Return the sampled dendrogram ;
(20)end if
5.2. Noisy Probability Calculation Method: PCBD

The noisy connection probability computing algorithm PCBD is shown in Algorithm 2. To avoid connection probability to be 0 or 1, we model as a distribution and prevent its expected value from becoming 0 or 1. When the Beta distribution has hyperparameters as α = β = 1, it corresponds to a uniform distribution. So, we treat the uniform distribution with the parameter as the Beta distribution with hyperparameters α = β = 1. In addition, the Beta distribution is conjugate with the Binomial distribution. When the conjugate prior distribution is the Beta distribution, the posteriori distribution is still the Beta distribution. And when the conjugate prior distribution is the uniform distribution with the parameter , the posteriori distribution is the Beta distribution as well. We could get the likelihood through Bayesian theory:

Input: evolving graphs , sampled dendrogram , privacy parameter , the length of the sliding window , ’s internal node r
Output: noisy probabilities ,
(1);
(2)compute the of and of ;
(3)for do
(4),;
(5), ;
(6)for end
(7),;
(8)
(9)for each r’child do
(10)  ;
(11)  PCBD ;
(12)end for

Given the evolving graphs in the sliding window , where is the length of the sliding window, we could compute the posteriori distribution by updating the hyperparameters. represents the maximum value of using Maximum A Posteriori (MAP). We update the hyperparameters as follows:

Thus, we could obtain from . To satisfy the differential privacy, we add noise to the hyperparameters and with application of the Laplace mechanism. The global sensitivity of PCBD is 2N, where N is the size of the node set which contains the nodes of all the graphs in the sliding window. The reason is that when we insert or delete any edge, the maximum effect it leads to is N. The sanitized hyperparameters are and , and

The PCBD algorithm is shown in Algorithm 2. In particular, we first get the union set which contains all the graphs in the sliding window and compute from . We should also compute from (lines 1-2). Then, we compute the increment of the hyperparameters and using and of each graph in the sliding window (lines 3-4). To satisfy the differential privacy, we add noise to and , and we get and (line 5). Next, we update the hyperparameters and and compute (lines 7-8). Finally, we use recursion method to compute the noisy probabilities of internal nodes of all the dendrograms.

5.3. Generation of Synthetic Graph

After getting the sampled dendrogram and the noisy probabilities , for each pair of nodes , we could find the common ancestor r of i, j which is closest to them. Then, we place an edge with independent probability and get the synthetic graph . Algorithm 3 shows the specific process.

Input: input graph , sampled dendrogram , privacy parameter
Output: synthetic graph
(1) root node of
(2) PCBD ;
(3)for each pair of vertices do
(4)  find the lowest common ancestor r of i and j in ;
(5)  place an edge in between i and j with independent probability
(6)for end
(7)return

6. Privacy Analysis of GHRG-DP

6.1. Sensitivity Computation of GHRG-DP

We first analyze the global sensitivity of DCBA algorithm. is the maximum change in the log-likelihood when any edge of graph misses. From [15], we know , where ,.

Theorem 4. and it is .

Proof. We first fix and treat as variable. Let , and we have . For all , we have . Then, we could know monotonically decreases. Since , we note that when is in (0,0.5), and when is in (0.5,1]. As has symmetric property, we have . When or 0, that is, or 1, gets the maximum value. Let , and we havewhere the constant . Let , we have , and .

6.2. Privacy Proof of GHRG-DP

Based on the sequential composition property, we prove that GHRG-DP satisfies -differential privacy.

Theorem 5. GHRG-DP satisfies -differential privacy.

Proof. We suppose that the length of evolving graphs is T, and there are m change points. Based on the sequential composition property, DCBA algorithm satisfies -differential privacy and PCBD algorithm satisfies -differential privacy. In Algorithm 3, we know that the synthetic graph is generated through sampled dendrogram and noisy probabilities, and it does not consume any privacy budget. Hence, GHRG-DP satisfies -differential privacy, and .

7. Experimental Results

7.1. Experiment Settings

We make use of two real-world evolving network datasets in our experiments: AS-caida [28] and AS-733. The dataset contains 122 CAIDA AS graphs, from January 2004 to November 2007. Each file contains a full AS graph derived from a set of RouteViews BGP table snapshots. In our experiment, we extracted snapshots weekly from Jan. 2004 to Apr. 2006. AS-733 collects the communication record of an Autonomous System (AS). It is collected from University of Oregon RouteViews Project. In our experiment, we extracted snapshots daily from Nov. 8, 1997, to Dec. 17, 1997. The statistics of portion data are shown in Tables 1 and 2.

We compare the accuracy of GHRG-DP with DDPA [29] and the straightforward method through (1) the relative error of degree distribution and shortest path length distribution (2) the running time. The relative error is represented as

As differential privacy needs to produce random noise, we measure the accuracy of the result by the median relative error where we run the Laplace mechanism 10 times.

7.2. Evaluation of GHRG-DP
7.2.1. Comparing with Other Approaches

To show the utility of GHRG-DP algorithm, we first compare the relative error of degree distribution and shortest path length distribution with DDPA and the straightforward method which we have mentioned above. From Figures 2 and 3, we can see that GHRG-DP outperforms the DDPA and straightforward method. This is because we add noise to the posterior parameters of the model instead of adding noise to the connection probabilities. GHRG-DP avoids perturbing the connection probabilities directly. Furthermore, we observe that, comparing the variance of the relative error during the entire time series, GHRG-DP also outperforms the straightforward method. This means that the output of GHRG-DP is more stable. This is because, during the process of PCBD algorithm, neighbor networks are similar to each other and it leads to the similar posterior parameters. Thus, the networks GHRG-DP generates are more stable.

Then, we compare the running time with DDPA and the straightforward method. From Table 3, we can see that GHRG-DP uses less time than DDPA and the straightforward method takes the most time. This is because, in GHRG-DP, we generate the dendrogram by adjustment in most time and it omits the MCMC process. However, the median of GHRG-DP is much less than the average, and the median of DDPA and the straightforward method is roughly equal to the average. This is because, when the result of change-point detection is false, the adjustment process of dendrogram omits the MCMC computing process. The running time is much less than the straightforward method. However, when the result of change-point detection is true, GHRG-DP needs to sample the dendrogram by MCMC. At this time, the running time is similar. As a result, the running time of GHRG-DP increases, and some part of running time becomes protruding.

7.2.2. Evaluation of GHRG-DP on Different

In Figures 4 and 5, we show how the assignment of privacy parameter affects the performance of GHRG-DP. We fix and assign for DCBA and PCBD, respectively. We compare the relative error of degree distribution and shortest path length distribution. From Figures 4 and 5, we can see, when is relatively large, the relative error always stays high. This is because, in PCBD, we perturb the model parameters and it affects the synthetic graph directly. This means that, compared with the process of sampling dendrogram under differential privacy, the process of perturbing the model parameters affects the performance of GHRG-DP much more.

7.2.3. Evaluation of GHRG-DP on Different Lengths of Sliding Window

In Figures 6 and 7, we show how the length of the sliding window affects the performance of GHRG-DP. We set the length to be , respectively. We compare the relative error of degree distribution and shortest path length distribution.

Form Figures 6 and 7, we can see, when the length of the sliding window is relatively long, the relative error always stays high. This is because, when we update the hyperparameters in PCBD, we need to compute the summation of and of all the nodes in the sliding window. It generates much cumulate noise.

8. Conclusions

In this paper, we investigate the problem of evolving network differential privacy. We first introduce a straightforward approach for publishing differentially private evolving network. This approach sanitizes each timestamp network, respectively. We observe that the amount of noise required in straightforward approach is really high. Thus, we propose our GHRG-DP algorithm. The main idea is adjusting the structure of dendrogram to reduce the amount of noise during the process of MCMC. We change the real value of connection to the distribution of connection probability. It satisfies the differential privacy through adding noise during the iteration process of updating the hyperparameters. It could fit the dynamic environment of evolving network. Formal privacy analysis and the results of extensive experiments show GHRG-DP can achieve a private evolving network with high utility. Our future research directions are the more effective change-point detection method in Algorithm 1 under differential privacy.

Data Availability

The data used to support the findings of this study can be accessed from http://snap.stanford.edu/.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This work was supported in part by the National Natural Science Foundation of China under grant numbers 61672179, 61370083, and 61402126, in part by the Natural Science Foundation of Heilongjiang Province under grant number F2015030, in part by the Science Foundation for Youths of Heilongjiang under grant number QC2016083, and in part by the Postdoctoral Foundation of Heilongjiang Province under grant number LBH-Z14071.